Introduction

Patterned distributions of individual plant species and species groups are likely to reflect patterns in environmental conditions and constraints imposed on plant survival, growth, and reproductive ability by abiotic and biotic factors, including physiological tolerance of levels of various environmental conditions (e.g., soil moisture, [Ca], wind speed), competition with other plants, and interactions with herbivores and pathogens. If variation in the abundance of focal species characterized by different growth forms or physiologies is correlated with environmental conditions, and not due merely to limited plant dispersal, it can suggest that those forms or physiologies are adapted to different states or levels of those conditions.

Sarracenia is a genus of 11 species of carnivorous plants, occurring on the coastal plains of the eastern United States, that produce pitcher-shaped leaves filled with liquid that function as pitfall traps.

Givnish (personal communication) observed a patchy distribution of S. alata in the extensive wet, sandy pine savannas at Buttercup Flats near Wiggins, Mississippi in the DeSoto National Forest. Such wet pine savannas of the Gulf Coastal plain of the United States contain an extraordinary diversity and abundance of carnivorous plants occurring down slope from mesic longleaf pine (Pinus palustris Mill.) dominated uplands (Brewer 1999). On these sites, S. alata tends to co-occur with crayfish mounds made of clay material (Givnish, pers. obs.), suggesting that the plant’s locally patchy distribution may have originated from a hidden variability in the physical composition of the substrate (Brewer 1999). Specifically, the presence of clay lenses in the subsoil may create localized anoxic zones that prevent nitrogen fixation by soil microbes, leading to wet, N-poor microsites that favor the establishment of carnivorous plants, due to low levels of soil nitrate and hostile, anoxic conditions that can work against root uptake of soil nutrients generally. I therefore aim to test the hypothesis that scattered but hidden clay lenses help generate the patchy distribution of Sarracenia, creating zones that are poor in nitrate and dissolved oxygen.

pattern1

pattern1

Thirteen transects crossing the described vegetation pattern were set up randomly within in the savanna at Buttercup Flats facing a compass bearing of 160 degrees. Each transect contained a pitcher-absent microsite in the middle and 10-m of pitcher-present microsites on each end. Transects 11, 12, and 13 intersected multiple pitcher-absent patches and were therefore evaluated at more locations compared to the other transects. Soil and tissue nutrients, ground Coverage of all species, and soil clay content clay content were measured at these sites using a stratified random design.

Analysis of Soil Nutrients

Plant Root Simulators (PRS, available from Western Ag) with cation and anion exchange resins be planted in pairs of three in two locations along the transect: 1) a random location at the pitcher-present microsite in the beginning ten meters of the transect and 2) a random location in each pitcher-absent microsite contained within the transect. Probe burial time was 46-47 days. All soil nutrients were measured in μg/cm2.

Each nutrient is evaluated independently for its association with microsite type (pitcher plant-present, pitcher plant-absent).

soilNPK = read.csv(file = "soilNPK.csv", header = TRUE)
soilNPK$Transect = factor(soilNPK$Transect) # interpreted as factor, not numerical values
str(soilNPK)
'data.frame':   31 obs. of  6 variables:
 $ Microsite: Factor w/ 2 levels "-","+": 2 1 2 1 2 1 2 1 2 1 ...
 $ Transect : Factor w/ 13 levels "1","2","3","4",..: 1 1 2 2 3 3 4 4 5 5 ...
 $ NO3      : num  0.6 5.2 0.9 0.9 3.8 2.8 3.2 1.5 1 1.7 ...
 $ NH4      : num  2.3 18.9 1.6 2.6 3.1 1.9 7.3 2 2.8 2.3 ...
 $ K        : num  128 321 105 130 147 ...
 $ P        : num  0.9 0.3 0.2 0.3 0.3 0.3 0.4 0.5 0.3 0.9 ...

Soil Nitrogren

Soil nitrogen is evaluated as both Total Nitrogen and its components:
Total Nitrogen is treated as the sum of soil nitrate (NO3-) and soil ammonium (NH4+).

Increased clay content would reduce drainage and create local flooding. The resulting anoxic environment would severely limit nitrification and favor denitrification, necesitating alternative nutrient acquisition strategies. We would, therefore, expect lower concentrations of soil nitrogen at clay-rich and (presumably) pitcher-present microsites (Craft and Chiang 2002)

Differences in soil levels of nitrate and ammonium between or within microsites might indicate a preferential aquisition of a given form of nitrogen based on physiological differences between dominant plants (Olde Venterink et al. 2003).

soilNPK$TN <- soilNPK$NO3 + soilNPK$NH4

Total Nitrogen

data visualization

layout(matrix(1:2,1,2))
hist(soilNPK$TN, main="", col="tan", breaks=20)
qqnorm(soilNPK$TN, main = '', ylab="soil K")
qqline(soilNPK$TN)

The histogram for total soil nitrogen concentrations is moderately right-skewed. Q-Q plots were used to test for normality because the Shapiro-Wilk test is not powerful at small sample size (n=31). Sample quantiles show noticible deviation from theoretical quantiles; Q-Q plot does suggests a significant deviation from normality. Two outliers (large TN values) can be seen in both the histogram and the Q-Q plot. they are on 2 different transects and represent different microsites.

subset(soilNPK, TN>15 )
   Microsite Transect NO3  NH4      K   P   TN
2          -        1 5.2 18.9 321.01 0.3 24.1
29         +       13 0.0 18.5 101.58 1.4 18.5

Let’s now just visualize the association between total soil nitrogen and microsite, and also see the pairing by transect (transect 1 is in red, transect 13 is in blue).

transect.color = rep("black",13)
transect.color[c(1,13)] = c("red", "blue")
layout(matrix(1:2,1,2))
with(soilNPK, interaction.plot(Microsite, Transect, TN, legend=FALSE, col=transect.color))
plot(TN~Microsite, data=soilNPK)

These two transects show large changes in mean TN (in opposite directions) between microsites (note that Transect 13 has several data points for the “-” microsite).

The outliers are much “less” outliers and the Q-Q plot shows no significant deviation from normality when we consider the log(potassium) values:

layout(matrix(1:4,2,2, byrow=T))
with(soilNPK, hist(log(TN), main="", col="tan", breaks=20))
with(soilNPK, qqnorm(log(TN), main = '', ylab="soil log(TN)"))
with(soilNPK, interaction.plot(Microsite, Transect, log(TN), legend=FALSE, col=transect.color))
plot(log(TN)~Microsite, data=soilNPK)

From this, TN values between microsites increase and decrease in an equal number of transects.

statistical analysis

If we analyze the log values, it seems that (at first sight) we can use a t-test or a 2-way ANOVA-based method. We will double-check the normality of residuals after we get these residuals.

At transects 11, 12, and 13, multiple data points (for soil nutrients only) were collected for pitcher plant-absent microsites. It is necessary to average these pitcher plant-absent data points for their respective transects to perform a t-test (which requires a single value per transect per microsite). This averaging is not necessary with a 2-way ANOVA-based approach.

A new vector is created specifically for the paired t-test.

table(soilNPK[1:2])
         Transect
Microsite 1 2 3 4 5 6 7 8 9 10 11 12 13
        - 1 1 1 1 1 1 1 1 1  1  4  2  2
        + 1 1 1 1 1 1 1 1 1  1  1  1  1
logTNmeans = with(soilNPK, tapply(log(TN), Microsite:Transect, mean))
TN_minus = logTNmeans[1:13]
TN_plus = logTNmeans[14:26]
TN_minus  # log values, in fact
      -:1       -:2       -:3       -:4       -:5       -:6       -:7 
3.1822118 1.2527630 1.5475625 1.2527630 1.3862944 2.0668628 1.0647107 
      -:8       -:9      -:10      -:11      -:12      -:13 
1.3350011 1.6486586 0.8754687 1.7639335 1.2668484 1.6134220 
TN_plus
      +:1       +:2       +:3       +:4       +:5       +:6       +:7 
1.0647107 0.9162907 1.9315214 2.3513753 1.3350011 1.5475625 1.7749524 
      +:8       +:9      +:10      +:11      +:12      +:13 
2.5952547 0.7884574 0.2623643 2.1860513 2.6390573 2.9177707 

A paired t-test is performed on pitcher plant-absent and pitcher plant-present TN soil concentrations.

t.test(TN_plus, TN_minus, paired = TRUE)

    Paired t-test

data:  TN_plus and TN_minus
t = 0.55005, df = 12, p-value = 0.5924
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.4678313  0.7838112
sample estimates:
mean of the differences 
              0.1579899 

TN concentrations are not significantly different between pitcher-present and pitcher-absent microsites (p=0.59).
We can do this again with the non-transformed TN values, but there is still no evidence of association (p=0.57).

TNmeans = with(soilNPK, tapply(TN, Microsite:Transect, mean))
TN_minus = TNmeans[1:13]
TN_plus = TNmeans[14:26]
t.test(TN_plus, TN_minus, paired=TRUE)

    Paired t-test

data:  TN_plus and TN_minus
t = 0.58788, df = 12, p-value = 0.5675
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -3.809498  6.624883
sample estimates:
mean of the differences 
               1.407692 

A linear model is also constructed as both a supplement and a comparison to the paired t-test. We get practically the same p-values, but we also get to see diagnostic plots to confirm that assumptions are met reasonably well.

The linear model is constructed with TN (or log(TN)) as the dependent variable, transect as the first factor, and microsite as the second factor. We first use the raw TN values, not transformed. The residuals look fairly normally distributed, but the first residual plot does not look great.
No evidence for association between TN and microsite (p=0.37).

fit <- lm(TN ~ Transect + Microsite, data=soilNPK)
anova(fit)
Analysis of Variance Table

Response: TN
          Df Sum Sq Mean Sq F value Pr(>F)
Transect  12 244.56  20.380  0.6534 0.7708
Microsite  1  26.81  26.813  0.8596 0.3668
Residuals 17 530.25  31.191               
drop1(fit, test="F")
Single term deletions

Model:
TN ~ Transect + Microsite
          Df Sum of Sq    RSS    AIC F value Pr(>F)
<none>                 530.25 116.02               
Transect  12   254.722 784.98 104.18  0.6805 0.7482
Microsite  1    26.813 557.07 115.55  0.8596 0.3668
layout(matrix(1:4,1,4, byrow=T))
plot(fit)

We can re-do this analysis, using log(TN) values this time. Still no evidence of an association (p=0.39). The residual plots look better (more homogeneous variance), so this analysis is preferred over the analysis of TN (untransformed).

fit <- lm(log(TN) ~ Transect + Microsite, data=soilNPK)
anova(fit)
Analysis of Variance Table

Response: log(TN)
          Df Sum Sq Mean Sq F value Pr(>F)
Transect  12 5.0649 0.42207  0.9186 0.5498
Microsite  1 0.3587 0.35873  0.7807 0.3892
Residuals 17 7.8111 0.45948               
drop1(fit, test="F")
Single term deletions

Model:
log(TN) ~ Transect + Microsite
          Df Sum of Sq     RSS     AIC F value Pr(>F)
<none>                  7.8111 -14.732               
Transect  12    5.2823 13.0934 -22.718  0.9580 0.5191
Microsite  1    0.3587  8.1699 -15.340  0.7807 0.3892
plot(fit)

Ammonium

data visualization

We repeat the same analysis using ammonium (NH4+) the cationic componenet of Total Nitrogen.

layout(matrix(1:2,1,2))
hist(soilNPK$NH4, main="", col="tan", breaks=20)
qqnorm(soilNPK$NH4, main = '', ylab="soil NH4")
qqline(soilNPK$NH4)

The histogran for soil ammonium concentrations is noticibly right-skewed. Q-Q plots were used to test for normality because the Shapiro-Wilk test is not powerful at small sample size (n=31). Sample quantiles do not noticible deviation from theoretical quantiles; the Q-Q plot suggests a significant deviation from normality. Four outliers (large ammonium values) can be seen in the Q-Q plots. they are on 4 different transects, representing both microsites equally.

subset(soilNPK, NH4>10 )
   Microsite Transect NO3  NH4      K   P   TN
2          -        1 5.2 18.9 321.01 0.3 24.1
25         -       11 0.0 12.2 110.25 0.8 12.2
26         +       12 0.0 14.0 111.88 0.8 14.0
29         +       13 0.0 18.5 101.58 1.4 18.5

Let’s now just visualize the association between ammonium and microsite, and also see the pairing by transect (transect 1 is in red, transects 12 and 13 are in blue).

transect.color = rep("black",13)
transect.color[c(1,13, 12)] = c("red", "blue", "blue")
layout(matrix(1:2,1,2))
with(soilNPK, interaction.plot(Microsite, Transect, NH4, legend=FALSE, col=transect.color))
plot(NH4~Microsite, data=soilNPK)

All transects, with the exeptions of the highlighted transects 1, 12, and 13, do not show large differences in ammonium. Ammonium concentrations tend to be slightly greater on the pitcher-present microsite of each transect. Q-Q plots show much less deviation from normality when we consider the log(ammonium) values:

layout(matrix(1:4,2,2, byrow=T))
with(soilNPK, hist(log(NH4), main="", col="tan", breaks=20))
with(soilNPK, qqnorm(log(NH4), main = '', ylab="soil log(NH4)"))
with(soilNPK, interaction.plot(Microsite, Transect, log(NH4), legend=FALSE, col=transect.color))
plot(log(NH4)~Microsite, data=soilNPK)

statistical analysis

If we analyze the log values, it seems that (at first sight) we can use a t-test or a 2-way ANOVA-based method. We will double-check the normality of residuals after we get these residuals.

At transects 11, 12, and 13, multiple data points (for soil nutrients only) were collected for pitcher plant-absent microsites. It is necessary to average these pitcher plant-absent data points for their respective transects to perform a t-test (which requires a single value per transect per microsite). This averaging is not necessary with a 2-way ANOVA-based approach.

A new vector is created specifically for the paired t-test.

table(soilNPK[1:2])
         Transect
Microsite 1 2 3 4 5 6 7 8 9 10 11 12 13
        - 1 1 1 1 1 1 1 1 1  1  4  2  2
        + 1 1 1 1 1 1 1 1 1  1  1  1  1
logNH4means = with(soilNPK, tapply(log(NH4), Microsite:Transect, mean))
NH4_minus = logNH4means[1:13]
NH4_plus = logNH4means[14:26]
NH4_minus  # log values, in fact
      -:1       -:2       -:3       -:4       -:5       -:6       -:7 
2.9391619 0.9555114 0.6418539 0.6931472 0.8329091 1.3609766 0.4054651 
      -:8       -:9      -:10      -:11      -:12      -:13 
0.2623643 0.9555114 0.3364722 1.6098012 0.6404669 1.2644618 
NH4_plus
      +:1       +:2       +:3       +:4       +:5       +:6       +:7 
0.8329091 0.4700036 1.1314021 1.9878743 1.0296194 1.3862944 0.6931472 
      +:8       +:9      +:10      +:11      +:12      +:13 
1.5040774 0.5306283 0.2623643 1.8082888 2.6390573 2.9177707 

A paired t-test is performed on pitcher plant-absent and pitcher plant-present K soil concentrations.

t.test(NH4_plus, NH4_minus, paired = TRUE)

    Paired t-test

data:  NH4_plus and NH4_minus
t = 1.1115, df = 12, p-value = 0.2881
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.3172923  0.9781129
sample estimates:
mean of the differences 
              0.3304103 

Ammonium concentrations are not significantly different between pitcher-present and pitcher-absent microsites (p=0.29).
We can do this again with the non-transformed Ammonium values, but there is still no evidence of association (p=0.48).

NH4means = with(soilNPK, tapply(NH4, Microsite:Transect, mean))
NH4_minus = NH4means[1:13]
NH4_plus = NH4means[14:26]
t.test(NH4_plus, NH4_minus, paired=TRUE)

    Paired t-test

data:  NH4_plus and NH4_minus
t = 0.72334, df = 12, p-value = 0.4833
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -2.991122  5.964199
sample estimates:
mean of the differences 
               1.486538 

A linear model is also constructed as both a supplement and a comparison to the paired t-test. We get practically the same p-values, but we also get to see diagnostic plots to confirm that assumptions are met reasonably well.

The linear model is constructed with ammonium (or log(ammonium)) as the dependent variable, transect as the first factor, and microsite as the second factor. We first use the raw ammonium values, not transformed. The residuals look fairly normally distributed, but the first residual plot and the Q-Q plot do not look great.
No evidence for association between ammonium and microsite (p=0.31).

fit <- lm(NH4 ~ Transect + Microsite, data=soilNPK)
anova(fit)
Analysis of Variance Table

Response: NH4
          Df Sum Sq Mean Sq F value Pr(>F)
Transect  12 219.84  18.320  0.7188 0.7159
Microsite  1  28.24  28.241  1.1081 0.3072
Residuals 17 433.27  25.487               
drop1(fit, test="F")
Single term deletions

Model:
NH4 ~ Transect + Microsite
          Df Sum of Sq    RSS     AIC F value Pr(>F)
<none>                 433.27 109.759               
Transect  12   236.380 669.65  99.256  0.7729 0.6700
Microsite  1    28.241 461.52 109.716  1.1081 0.3072
layout(matrix(1:4,1,4, byrow=T))
plot(fit)

We can re-do this analysis, using log(ammonium) values this time. Still no evidence of an association (p=0.17). The residual plots look better (more homogeneous variance) and the Q-Q plot show less deviation from normality, so this analysis is preferred over the analysis of ammonium (untransformed).

fit <- lm(log(NH4) ~ Transect + Microsite, data=soilNPK)
anova(fit)
Analysis of Variance Table

Response: log(NH4)
          Df Sum Sq Mean Sq F value Pr(>F)
Transect  12 7.1205 0.59338  1.1050 0.4147
Microsite  1 1.1116 1.11155  2.0699 0.1684
Residuals 17 9.1290 0.53700               
drop1(fit, test="F")
Single term deletions

Model:
log(NH4) ~ Transect + Microsite
          Df Sum of Sq    RSS      AIC F value Pr(>F)
<none>                  9.129  -9.8983               
Transect  12    7.8262 16.955 -14.7057  1.2145 0.3480
Microsite  1    1.1116 10.241  -8.3364  2.0699 0.1684
plot(fit)

Nitrate

We, again, repeat the same analysis, this time using nitrate (NO3-) the anionic componenet of Total Nitrogen.

data visualization

layout(matrix(1:2,1,2))
hist(soilNPK$NO3, main="", col="tan", breaks=20)
qqnorm(soilNPK$NO3, main = '', ylab="soil NO3")
qqline(soilNPK$NO3)

The histogran for soil nitrate concentrations is moderately right-skewed. Q-Q plots were used to test for normality because the Shapiro-Wilk test is not powerful at small sample size (n=31). Sample quantiles do not show noticible deviation from Q-Q plot does not suggest any significant deviations from normality (though the first few quantiles look a little weird). Two outliers (large nitrate values) can be seen in the Q-Q plots. they are on 2 different transects and represent different microsites.

subset(soilNPK, NO3>4 )
   Microsite Transect NO3  NH4      K   P   TN
2          -        1 5.2 18.9 321.01 0.3 24.1
15         +        8 8.9  4.5  83.52 0.7 13.4

Let’s now just visualize the association between nitrate and microsite, and also see the pairing by transect (transect 1 is in red, transect 8 is in blue).

transect.color = rep("black",13)
transect.color[c(1,8)] = c("red", "blue")
layout(matrix(1:2,1,2))
with(soilNPK, interaction.plot(Microsite, Transect, NO3, legend=FALSE, col=transect.color))
plot(NO3~Microsite, data=soilNPK)

Two large changes in nitrate concentrates (those of the highlighted transects) are visible, going in different directions, between microsites.

The outliers are much “less” outliers when we consider the squareroot(nitrate) values. A square root tranformation was chosen for nitrate due to the presence of several zeroes in the data. The Q-Q plot has lost its outliers but still does not look great:

layout(matrix(1:4,2,2, byrow=T))
with(soilNPK, hist(sqrt(NO3), main="", col="tan", breaks=20))
with(soilNPK, qqnorm(sqrt(NO3), main = '', ylab="soil log(NO3)"))
with(soilNPK, interaction.plot(Microsite, Transect, sqrt(NO3), legend=FALSE, col=transect.color))
plot(sqrt(NO3)~Microsite, data=soilNPK)

From this, most transects show slight changes in nitrate concentration across microsites. These changes go in different directions in an equal frequency.

statistical analysis

If we analyze the squart root values, it seems that (at first sight) we can use a t-test or a 2-way ANOVA-based method. We will double-check the normality of residuals after we get these residuals.

At transects 11, 12, and 13, multiple data points (for soil nutrients only) were collected for pitcher plant-absent microsites. It is necessary to average these pitcher plant-absent data points for their respective transects to perform a t-test (which requires a single value per transect per microsite). This averaging is not necessary with a 2-way ANOVA-based approach.

A new vector is created specifically for the paired t-test.

table(soilNPK[1:2])
         Transect
Microsite 1 2 3 4 5 6 7 8 9 10 11 12 13
        - 1 1 1 1 1 1 1 1 1  1  4  2  2
        + 1 1 1 1 1 1 1 1 1  1  1  1  1
sqrtNO3means = with(soilNPK, tapply(sqrt(NO3), Microsite:Transect, mean))
NO3_minus = sqrtNO3means[1:13]
NO3_plus = sqrtNO3means[14:26]
NO3_minus  # log values, in fact
      -:1       -:2       -:3       -:4       -:5       -:6       -:7 
2.2803509 0.9486833 1.6733201 1.2247449 1.3038405 2.0000000 1.1832160 
      -:8       -:9      -:10      -:11      -:12      -:13 
1.5811388 1.6124515 1.0000000 0.5049690 1.2831928 1.2159615 
NO3_plus
      +:1       +:2       +:3       +:4       +:5       +:6       +:7 
0.7745967 0.9486833 1.9493589 1.7888544 1.0000000 0.8366600 1.9748418 
      +:8       +:9      +:10      +:11      +:12      +:13 
2.9832868 0.7071068 0.0000000 1.6733201 0.0000000 0.0000000 

A paired t-test is performed on pitcher plant-absent and pitcher plant-present K soil concentrations.

t.test(NO3_plus, NO3_minus, paired = TRUE)

    Paired t-test

data:  NO3_plus and NO3_minus
t = -0.87215, df = 12, p-value = 0.4002
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.8544107  0.3659245
sample estimates:
mean of the differences 
             -0.2442431 

K concentrations are not significantly different between pitcher-present and pitcher-absent microsites (p=0.92).
We can do this again with the non-transformed K values, but there is still no evidence of association (p=0.92).

NO3means = with(soilNPK, tapply(NO3, Microsite:Transect, mean))
NO3_minus = NO3means[1:13]
NO3_plus = NO3means[14:26]
t.test(NO3_plus, NO3_minus, paired=TRUE)

    Paired t-test

data:  NO3_plus and NO3_minus
t = -0.09884, df = 12, p-value = 0.9229
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -1.816912  1.659220
sample estimates:
mean of the differences 
            -0.07884615 

A linear model is also constructed as both a supplement and a comparison to the paired t-test. We get practically the same p-values, but we also get to see diagnostic plots to confirm that assumptions are met reasonably well.

The linear model is constructed with nitrate (or squareroot(nitrte)) as the dependent variable, transect as the first factor, and microsite as the second factor. We first use the raw nitrate values, not transformed. The residuals look fairly normally distributed, and the Q-Q plot suggests normality.
No evidence for association between nitrate and microsite (p=0.28).

fit <- lm(NO3 ~ Transect + Microsite, data=soilNPK)
anova(fit)
Analysis of Variance Table

Response: NO3
          Df Sum Sq Mean Sq F value Pr(>F)
Transect  12 51.902  4.3252  1.3541 0.2765
Microsite  1  0.019  0.0185  0.0058 0.9402
Residuals 17 54.301  3.1942               
drop1(fit, test="F")
Single term deletions

Model:
NO3 ~ Transect + Microsite
          Df Sum of Sq     RSS    AIC F value Pr(>F)
<none>                  54.301 45.377               
Transect  12    51.487 105.788 42.051  1.3432 0.2815
Microsite  1     0.019  54.320 43.388  0.0058 0.9402
layout(matrix(1:4,1,4, byrow=T))
plot(fit)

We can re-do this analysis, using squareroot(nitrate) values this time. Still no evidence of an association (p=0.34). The residual plots look (slightly) better (more homogeneous variance), so this analysis is preferred over the analysis of nitrate (untransformed).

fit <- lm(sqrt(NO3) ~ Transect + Microsite, data=soilNPK)
anova(fit)
Analysis of Variance Table

Response: sqrt(NO3)
          Df Sum Sq Mean Sq F value Pr(>F)
Transect  12 6.8471 0.57059  1.1713 0.3732
Microsite  1 0.3833 0.38333  0.7869 0.3874
Residuals 17 8.2817 0.48716               
drop1(fit, test="F")
Single term deletions

Model:
sqrt(NO3) ~ Transect + Microsite
          Df Sum of Sq     RSS     AIC F value Pr(>F)
<none>                  8.2817 -12.918               
Transect  12    7.1737 15.4555 -17.577  1.2271 0.3409
Microsite  1    0.3833  8.6651 -13.515  0.7869 0.3874
plot(fit)

Soil Phosphorus

Soil phosphorus, in contrast to nitrogen, is expected to increase with local inundation, which chemically reduces Fe bound P, releasing Phosphate into the soil (Craft and Chiang 2002).

data visualization

layout(matrix(1:2,1,2))
hist(soilNPK$P, main="", col="tan", breaks=20)
qqnorm(soilNPK$P, main = '', ylab="soil P")
qqline(soilNPK$P)

The histogran for soil phosphorus concentrations is slightly left-skewed. Q-Q plots were used to test for normality because the Shapiro-Wilk test is not powerful at small sample size (n=31). Sample quantiles show noticible deviation from theoretical quantiles; the Q-Q plot suggests significant deviations from normality. One outlier (large P value) can be seen in the Q-Q plots. This data point is on transect 11 and in the “+” microsite: pitcher plants present

subset(soilNPK, P>1.4 )
   Microsite Transect NO3 NH4      K   P  TN
11         +        6 0.7   4 103.23 1.6 4.7

Let’s now just visualize the association between phosphorus and microsite, and also see the pairing by transect (transect 9 is in red, transects 6 and 13 is in blue).

transect.color = rep("black",13)
transect.color[c(9,13,6)] = c("red", "blue", "blue")
layout(matrix(1:2,1,2))
with(soilNPK, interaction.plot(Microsite, Transect, P, legend=FALSE, col=transect.color))
plot(P~Microsite, data=soilNPK)

The transects with the largest changes in soil (6,9,13) phosphorus do not correspond to the outlier (Transect 11). Soil phosphorus, in general, increases at the pitcher-present site.

The outliers are much “less” outliers when we consider the log(phosphorus) values (but the Q-Q plot still does not look great!):

layout(matrix(1:4,2,2, byrow=T))
with(soilNPK, hist(log(P), main="", col="tan", breaks=20))
with(soilNPK, qqnorm(log(P), main = '', ylab="soil log(P)"))
with(soilNPK, interaction.plot(Microsite, Transect, log(P), legend=FALSE, col=transect.color))
plot(log(P)~Microsite, data=soilNPK)

statistical analysis

If we analyze the log values, it seems that (at first sight) we can use a t-test or a 2-way ANOVA-based method. We will double-check the normality of residuals after we get these residuals.

At transects 11, 12, and 13, multiple data points (for soil nutrients only) were collected for pitcher plant-absent microsites. It is necessary to average these pitcher plant-absent data points for their respective transects to perform a t-test (which requires a single value per transect per microsite). This averaging is not necessary with a 2-way ANOVA-based approach.

A new vector is created specifically for the paired t-test.

table(soilNPK[1:2])
         Transect
Microsite 1 2 3 4 5 6 7 8 9 10 11 12 13
        - 1 1 1 1 1 1 1 1 1  1  4  2  2
        + 1 1 1 1 1 1 1 1 1  1  1  1  1
logPmeans = with(soilNPK, tapply(log(P), Microsite:Transect, mean))
P_minus = logPmeans[1:13]
P_plus = logPmeans[14:26]
P_minus  # log values, in fact
        -:1         -:2         -:3         -:4         -:5         -:6 
-1.20397280 -1.20397280 -1.20397280 -0.69314718 -0.10536052 -0.51082562 
        -:7         -:8         -:9        -:10        -:11        -:12 
 0.18232156 -1.20397280  0.18232156  0.09531018 -0.18600488  0.13118213 
       -:13 
-1.06013177 
P_plus
        +:1         +:2         +:3         +:4         +:5         +:6 
-0.10536052 -1.60943791 -1.20397280 -0.91629073 -1.20397280  0.47000363 
        +:7         +:8         +:9        +:10        +:11        +:12 
-0.22314355 -0.35667494 -0.69314718  0.09531018 -0.69314718 -0.22314355 
       +:13 
 0.33647224 

A paired t-test is performed on pitcher plant-absent and pitcher plant-present K soil concentrations.

t.test(P_plus, P_minus, paired = TRUE)

    Paired t-test

data:  P_plus and P_minus
t = 0.15858, df = 12, p-value = 0.8766
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.4446241  0.5144272
sample estimates:
mean of the differences 
             0.03490159 

P concentrations are not significantly different between pitcher-present and pitcher-absent microsites (p=0.88).
We can do this again with the non-transformed P values, but there is still no evidence of association (p=0.85).

Pmeans = with(soilNPK, tapply(P, Microsite:Transect, mean))
P_minus = Pmeans[1:13]
P_plus = Pmeans[14:26]
t.test(P_plus, P_minus, paired=TRUE)

    Paired t-test

data:  P_plus and P_minus
t = 0.1976, df = 12, p-value = 0.8467
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.3085107  0.3700491
sample estimates:
mean of the differences 
             0.03076923 

A linear model is also constructed as both a supplement and a comparison to the paired t-test. We get practically the same p-values, but we also get to see diagnostic plots to confirm that assumptions are met reasonably well.

The linear model is constructed with P (or log(P)) as the dependent variable, transect as the first factor, and microsite as the second factor. We first use the raw K values, not transformed. The residuals look fairly normally distributed.
No evidence for association between K and microsite (p=0.84).

fit <- lm(P ~ Transect + Microsite, data=soilNPK)
anova(fit)
Analysis of Variance Table

Response: P
          Df Sum Sq  Mean Sq F value Pr(>F)
Transect  12 2.1714 0.180951  1.2977 0.3036
Microsite  1 0.0058 0.005796  0.0416 0.8409
Residuals 17 2.3705 0.139443               
drop1(fit, test="F")
Single term deletions

Model:
P ~ Transect + Microsite
          Df Sum of Sq    RSS     AIC F value Pr(>F)
<none>                 2.3705 -51.697               
Transect  12    2.1772 4.5477 -55.500  1.3011 0.3019
Microsite  1    0.0058 2.3763 -53.621  0.0416 0.8409
layout(matrix(1:4,1,4, byrow=T))
plot(fit)

We can re-do this analysis, using log(P) values this time. Still no evidence of an association (p=0.86). The residual plots look better (more homogeneous variance), so this analysis is preferred over the analysis of P (untransformed).

fit <- lm(log(P) ~ Transect + Microsite, data=soilNPK)
anova(fit)
Analysis of Variance Table

Response: log(P)
          Df Sum Sq Mean Sq F value Pr(>F)
Transect  12 5.8384 0.48653  1.6440 0.1696
Microsite  1 0.0087 0.00865  0.0292 0.8663
Residuals 17 5.0310 0.29594               
drop1(fit, test="F")
Single term deletions

Model:
log(P) ~ Transect + Microsite
          Df Sum of Sq     RSS     AIC F value Pr(>F)
<none>                  5.0310 -28.369               
Transect  12    5.8414 10.8724 -28.481  1.6449 0.1694
Microsite  1    0.0087  5.0396 -30.316  0.0292 0.8663
plot(fit)

Soil Potassium

Soil Potassium may be less abundant at microsites with lower clay content due a probable increase in nutrient leeching (Inglett, Reddy, and Corstanje 2005).

data visualization

layout(matrix(1:2,1,2))
hist(soilNPK$K, main="", col="tan", breaks=20)
qqnorm(soilNPK$K, main = '', ylab="soil K")
qqline(soilNPK$K)

The histogran for soil potassium concentrations is slightly right-skewed. Q-Q plots were used to test for normality because the Shapiro-Wilk test is not powerful at small sample size (n=31). Sample quantiles do not show noticible deviation from theoretical quantiles; Q-Q plot does not suggest any significant deviations from normality. Three outliers (large K values) can be seen in the Q-Q plots. they are on 3 different transects, but all in the “-” microsite: pitcher plants absent (note that one of these transects, #13, has another observation for “-” microsite)

subset(soilNPK, K>200 )
   Microsite Transect NO3  NH4      K   P   TN
2          -        1 5.2 18.9 321.01 0.3 24.1
8          -        4 1.5  2.0 229.97 0.5  3.5
31         -       13 2.2  3.8 225.55 0.6  6.0

Let’s now just visualize the association between potassium and microsite, and also see the pairing by transect (transect 1 is in red, transect 4 is in blue).

transect.color = rep("black",13)
transect.color[c(1,4)] = c("red", "blue")
layout(matrix(1:2,1,2))
with(soilNPK, interaction.plot(Microsite, Transect, K, legend=FALSE, col=transect.color))
plot(K~Microsite, data=soilNPK)

The outliers are much “less” outliers when we consider the log(potassium) values:

layout(matrix(1:4,2,2, byrow=T))
with(soilNPK, hist(log(K), main="", col="tan", breaks=20))
with(soilNPK, qqnorm(log(K), main = '', ylab="soil log(K)"))
with(soilNPK, interaction.plot(Microsite, Transect, log(K), legend=FALSE, col=transect.color))
plot(log(K)~Microsite, data=soilNPK)

From this, 4 transects have K higher in the pitcher-present area, and 9 transects have K higher in the pitcher-absent area.

statistical analysis

If we analyze the log values, it seems that (at first sight) we can use a t-test or a 2-way ANOVA-based method. We will double-check the normality of residuals after we get these residuals.

At transects 11, 12, and 13, multiple data points (for soil nutrients only) were collected for pitcher plant-absent microsites. It is necessary to average these pitcher plant-absent data points for their respective transects to perform a t-test (which requires a single value per transect per microsite). This averaging is not necessary with a 2-way ANOVA-based approach.

A new vector is created specifically for the paired t-test.

table(soilNPK[1:2])
         Transect
Microsite 1 2 3 4 5 6 7 8 9 10 11 12 13
        - 1 1 1 1 1 1 1 1 1  1  4  2  2
        + 1 1 1 1 1 1 1 1 1  1  1  1  1
logKmeans = with(soilNPK, tapply(log(K), Microsite:Transect, mean))
K_minus = logKmeans[1:13]
K_plus = logKmeans[14:26]
K_minus  # log values, in fact
     -:1      -:2      -:3      -:4      -:5      -:6      -:7      -:8 
5.771472 4.868534 4.652435 5.437949 4.959482 3.810876 4.453417 4.704925 
     -:9     -:10     -:11     -:12     -:13 
4.279994 4.211387 4.529404 4.798584 4.838847 
K_plus
     +:1      +:2      +:3      +:4      +:5      +:6      +:7      +:8 
4.850858 4.654817 4.988935 4.259718 4.596129 4.636960 4.755915 4.425086 
     +:9     +:10     +:11     +:12     +:13 
4.076859 4.863913 4.371218 4.717427 4.620847 

A paired t-test is performed on pitcher plant-absent and pitcher plant-present K soil concentrations.

t.test(K_plus, K_minus, paired = TRUE)

    Paired t-test

data:  K_plus and K_minus
t = -0.74329, df = 12, p-value = 0.4716
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.4531952  0.2226375
sample estimates:
mean of the differences 
             -0.1152788 

K concentrations are not significantly different between pitcher-present and pitcher-absent microsites (p=0.47).
We can do this again with the non-transformed K values, but there is still no evidence of association (p=0.23).

Kmeans = with(soilNPK, tapply(K, Microsite:Transect, mean))
K_minus = Kmeans[1:13]
K_plus = Kmeans[14:26]
t.test(K_plus, K_minus, paired=TRUE)

    Paired t-test

data:  K_plus and K_minus
t = -1.2574, df = 12, p-value = 0.2325
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -72.50872  19.44179
sample estimates:
mean of the differences 
              -26.53346 

A linear model is also constructed as both a supplement and a comparison to the paired t-test. We get practically the same p-values, but we also get to see diagnostic plots to confirm that assumptions are met reasonably well.

The linear model is constructed with K (or log(K)) as the dependent variable, transect as the first factor, and microsite as the second factor. We first use the raw K values, not transformed. The residuals look fairly normally distributed, but the first residual plot does not look great.
No evidence for association between K and microsite (p=0.22).

fit <- lm(K ~ Transect + Microsite, data=soilNPK)
anova(fit)
Analysis of Variance Table

Response: K
          Df Sum Sq Mean Sq F value Pr(>F)
Transect  12  39671  3306.0  1.0959 0.4207
Microsite  1   4940  4939.6  1.6375 0.2179
Residuals 17  51282  3016.6               
drop1(fit, test="F")
Single term deletions

Model:
K ~ Transect + Microsite
          Df Sum of Sq   RSS    AIC F value Pr(>F)
<none>                 51282 257.75               
Transect  12     40938 92220 251.94  1.1309 0.3981
Microsite  1      4940 56222 258.60  1.6375 0.2179
layout(matrix(1:4,1,4, byrow=T))
plot(fit)

We can re-do this analysis, using log(K) values this time. Still no evidence of an association (p=0.47). The residual plots look better (more homogeneous variance), so this analysis is preferred over the analysis of K (untransformed).

fit <- lm(log(K) ~ Transect + Microsite, data=soilNPK)
anova(fit)
Analysis of Variance Table

Response: log(K)
          Df  Sum Sq Mean Sq F value Pr(>F)
Transect  12 2.12012 0.17668  0.9749 0.5064
Microsite  1 0.10048 0.10048  0.5544 0.4667
Residuals 17 3.08095 0.18123               
drop1(fit, test="F")
Single term deletions

Model:
log(K) ~ Transect + Microsite
          Df Sum of Sq    RSS     AIC F value Pr(>F)
<none>                 3.0810 -43.571               
Transect  12   2.15191 5.2329 -51.150  0.9895 0.4954
Microsite  1   0.10048 3.1814 -44.576  0.5544 0.4667
plot(fit)

Analysis of Sedge tissue nutrients

Aboveground biomass of 6 species of sedges in the Rhynchospora genus (R. megalocarpa, R. compressa, R. gracilenta, R. plumosa, R. rariflora, R. harveyi var harveyi) were collected, dried, and analyzed for % N, P, and K content by the University of Wisconsin Soil and Forage Lab, as a measure of the availability of these nutrients to non-carnivorous plants. All species are C3.

All tissue nutrients were measured in % of Aboveground Dry Mass.

tissuenutrients <- read.table("tissueNPK.csv", header = TRUE, sep=",")
tissuenutrients$Transect = factor(tissuenutrients$Transect)
str(tissuenutrients)
'data.frame':   29 obs. of  5 variables:
 $ Microsite: Factor w/ 2 levels "-","+": 2 1 2 1 2 1 2 1 2 1 ...
 $ Transect : Factor w/ 13 levels "1","2","3","4",..: 1 1 2 2 3 3 4 4 5 5 ...
 $ TN       : num  1.55 1.07 0.63 0.7 1.09 0.81 1.01 0.52 1.08 0.55 ...
 $ P        : num  0.14 0.0869 0.08 0.06 0.0868 ...
 $ K        : num  0.67 0.59 0.81 0.61 0.52 0.89 0.6 0.63 0.61 0.9 ...

Tissue Nitrogen

Tissue Nitrogen was analyzed only as Total Nitrogen.

Lower soil nitrogen concentrations, combined with inhibition of root uptake on anoxic microsites, could likely result in lower plant tissue N on clay-rich microsites (Olde Venterink et al. 2003).

data visualization

layout(matrix(1:2,1,2))
hist(tissuenutrients$TN, main="", col="tan", breaks=20)
qqnorm(tissuenutrients$TN, main = '', ylab="soil TN")
qqline(tissuenutrients$TN)

The histogram for tissue nitrogen concentrations is moderately right-skewed unimodal. Q-Q plots were used to test for normality because the Shapiro-Wilk test is not powerful at small sample size (n=29). Sample quantiles do not show noticible deviation from theoretical quantiles; the Q-Q plot does not suggest any significant deviations from normality. One outliers (a large TN value) can be seen in the histogram and Q-Q plot. This data point is on Transect 1 and in the “+” microsite: pitcher plants present

subset(tissuenutrients, TN>1.4 )
  Microsite Transect   TN    P    K
1         +        1 1.55 0.14 0.67

Let’s now just visualize the association between tissue nitrogen and microsite, and also see the pairing by transect (transect 1 is in red, transect 13 is in blue).

transect.color = rep("black",13)
transect.color[c(1,13)] = c("red", "blue")
layout(matrix(1:2,1,2))
with(tissuenutrients, interaction.plot(Microsite, Transect, TN, legend=FALSE, col=transect.color))
plot(TN~Microsite, data=tissuenutrients)

The outliers are much “less” outliers when we consider the log(potassium) values:

layout(matrix(1:4,2,2, byrow=T))
with(tissuenutrients, hist(log(TN), main="", col="tan", breaks=20))
with(tissuenutrients, qqnorm(log(TN), main = '', ylab="soil log(TN)"))
with(tissuenutrients, interaction.plot(Microsite, Transect, log(TN), legend=FALSE, col=transect.color))
plot(log(TN)~Microsite, data=tissuenutrients)

From this, TN concentrations tend to change between microsites roughly equally in both directons.

statistical analysis

If we analyze the log values, it seems that (at first sight) we can use a t-test or a 2-way ANOVA-based method. We will double-check the normality of residuals after we get these residuals.

At transects 11, 12, and 13, multiple data points (for soil nutrients only) were collected for pitcher plant-absent microsites. It is necessary to average these pitcher plant-absent data points for their respective transects to perform a t-test (which requires a single value per transect per microsite). This averaging is not necessary with a 2-way ANOVA-based approach.

A new vector is created specifically for the paired t-test.

table(tissuenutrients[1:2])
         Transect
Microsite 1 2 3 4 5 6 7 8 9 10 11 12 13
        - 1 1 1 1 1 1 1 1 1  1  4  1  1
        + 1 1 1 1 1 1 1 1 1  1  1  1  1
logTNmeans = with(tissuenutrients, tapply(log(TN), Microsite:Transect, mean))
TN_minus = logTNmeans[1:13]
TN_plus = logTNmeans[14:26]
TN_minus  # log values, in fact
        -:1         -:2         -:3         -:4         -:5         -:6 
 0.06765865 -0.35667494 -0.21072103 -0.65392647 -0.59783700  0.03922071 
        -:7         -:8         -:9        -:10        -:11        -:12 
 0.02955880 -0.06187540 -0.02020271 -0.22314355  0.04306491 -0.31471074 
       -:13 
 0.15700375 
TN_plus
         +:1          +:2          +:3          +:4          +:5 
 0.438254931 -0.462035460  0.086177696  0.009950331  0.076961041 
         +:6          +:7          +:8          +:9         +:10 
 0.076961041  0.009950331 -0.030459207  0.086177696  0.122217633 
        +:11         +:12         +:13 
-0.020202707 -0.385662481 -0.301105093 

A paired t-test is performed on pitcher plant-absent and pitcher plant-present TN tissue concentrations.

t.test(TN_plus, TN_minus, paired = TRUE)

    Paired t-test

data:  TN_plus and TN_minus
t = 1.5659, df = 12, p-value = 0.1434
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.05449268  0.33291895
sample estimates:
mean of the differences 
              0.1392131 

TN tissue concentrations are not significantly different between pitcher-present and pitcher-absent microsites (p=0.14).
We can do this again with the non-transformed TN tissue values, but there is still no evidence of association (p=0.12).

TNmeans = with(tissuenutrients, tapply(TN, Microsite:Transect, mean))
TN_minus = TNmeans[1:13]
TN_plus = TNmeans[14:26]
t.test(TN_plus, TN_minus, paired=TRUE)

    Paired t-test

data:  TN_plus and TN_minus
t = 1.6338, df = 12, p-value = 0.1282
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.04234009  0.29618624
sample estimates:
mean of the differences 
              0.1269231 

A linear model is also constructed as both a supplement and a comparison to the paired t-test. We get practically the same p-values, but we also get to see diagnostic plots to confirm that assumptions are met reasonably well.

The linear model is constructed with TN (or log(TN)) as the dependent variable, transect as the first factor, and microsite as the second factor. We first use the raw TN values, not transformed. The residuals look fairly normally distributed, but the first residual plot does not look great.
No evidence for association between TN and microsite (p=0.12).

fit <- lm(TN ~ Transect + Microsite, data=tissuenutrients)
anova(fit)
Analysis of Variance Table

Response: TN
          Df  Sum Sq  Mean Sq F value Pr(>F)
Transect  12 0.73255 0.061046  1.7249 0.1583
Microsite  1 0.09506 0.095061  2.6861 0.1220
Residuals 15 0.53086 0.035391               
drop1(fit, test="F")
Single term deletions

Model:
TN ~ Transect + Microsite
          Df Sum of Sq     RSS     AIC F value Pr(>F)
<none>                 0.53086 -88.016               
Transect  12   0.76379 1.29465 -86.163  1.7985 0.1409
Microsite  1   0.09506 0.62592 -85.239  2.6861 0.1220
layout(matrix(1:4,1,4, byrow=T))
plot(fit)

We can re-do this analysis, using log(TN) values this time. Still no evidence of an association (p=0.13). The residual plots look better (more homogeneous variance), so this analysis is preferred over the analysis of TN (untransformed).

fit <- lm(log(TN) ~ Transect + Microsite, data=tissuenutrients)
anova(fit)
Analysis of Variance Table

Response: log(TN)
          Df  Sum Sq  Mean Sq F value Pr(>F)
Transect  12 0.91612 0.076343  1.6971 0.1655
Microsite  1 0.11542 0.115416  2.5657 0.1301
Residuals 15 0.67477 0.044984               
drop1(fit, test="F")
Single term deletions

Model:
log(TN) ~ Transect + Microsite
          Df Sum of Sq     RSS     AIC F value Pr(>F)
<none>                 0.67477 -81.060               
Transect  12   0.95864 1.63341 -79.422  1.7759 0.1460
Microsite  1   0.11542 0.79018 -78.481  2.5657 0.1301
plot(fit)

Tissue Phosphorus

Increased soil P avaliablilty due to frequent inundation would likely result in increased tissue P on clay-rich microsites (Craft and Chiang 2002)

data visualization

layout(matrix(1:2,1,2))
hist(tissuenutrients$P, main="", col="tan", breaks=20)
qqnorm(tissuenutrients$P, main = '', ylab="soil P")
qqline(tissuenutrients$P)

The histogram for tissue phosphorus concentrations is slightly right-skewed. Q-Q plots were used to test for normality because the Shapiro-Wilk test is not powerful at small sample size (n=31). Sample quantiles do not show noticible deviation from Q-Q plot does not suggest any significant deviations from normality. One outlier (large P value) can be seen in the Q-Q plot and histogram. This data point is on Transect 1 in the “+” microsite: pitcher plants absent

subset(tissuenutrients, P>0.12 )
  Microsite Transect   TN    P    K
1         +        1 1.55 0.14 0.67

Let’s now just visualize the association between potassium and microsite, and also see the pairing by transect (transect 1 is in red, transects 6 an 13 are in blue).

transect.color = rep("black",13)
transect.color[c(1,13,6)] = c("red", "blue", "blue")
layout(matrix(1:2,1,2))
with(tissuenutrients, interaction.plot(Microsite, Transect, P, legend=FALSE, col=transect.color))
plot(K~Microsite, data=tissuenutrients)

The outliers are much “less” outliers when we consider the log(phosphorus) values:

layout(matrix(1:4,2,2, byrow=T))
with(tissuenutrients, hist(log(P), main="", col="tan", breaks=20))
with(tissuenutrients, qqnorm(log(P), main = '', ylab="soil log(P)"))
with(tissuenutrients, interaction.plot(Microsite, Transect, log(P), legend=FALSE, col=transect.color))
plot(log(P)~Microsite, data=tissuenutrients)

From this, tissue P concentrations generally decrease in the pitcher-present microsite.

statistical analysis

If we analyze the log values, it seems that (at first sight) we can use a t-test or a 2-way ANOVA-based method. We will double-check the normality of residuals after we get these residuals.

table(tissuenutrients[1:2])
         Transect
Microsite 1 2 3 4 5 6 7 8 9 10 11 12 13
        - 1 1 1 1 1 1 1 1 1  1  4  1  1
        + 1 1 1 1 1 1 1 1 1  1  1  1  1
logPmeans = with(tissuenutrients, tapply(log(P), Microsite:Transect, mean))
P_minus = logPmeans[1:13]
P_plus = logPmeans[14:26]
P_minus  # log values, in fact
      -:1       -:2       -:3       -:4       -:5       -:6       -:7 
-2.442560 -2.813411 -2.323911 -2.302585 -2.318024 -2.207275 -2.578154 
      -:8       -:9      -:10      -:11      -:12      -:13 
-2.579459 -2.575116 -2.570648 -2.574808 -2.302585 -2.317587 
P_plus
      +:1       +:2       +:3       +:4       +:5       +:6       +:7 
-1.966113 -2.525729 -2.444633 -2.315498 -2.562315 -2.564184 -2.441365 
      +:8       +:9      +:10      +:11      +:12      +:13 
-2.449069 -2.445175 -2.574788 -2.569667 -2.659260 -2.723707 

A paired t-test is performed on pitcher plant-absent and pitcher plant-present P tissue concentrations.

t.test(P_plus, P_minus, paired = TRUE)

    Paired t-test

data:  P_plus and P_minus
t = -0.34997, df = 12, p-value = 0.7324
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.1864111  0.1348144
sample estimates:
mean of the differences 
            -0.02579836 

P concentrations are not significantly different between pitcher-present and pitcher-absent microsites (p=0.73).
We can do this again with the non-transformed K values, but there is still no evidence of association (p=0.78).

Pmeans = with(tissuenutrients, tapply(P, Microsite:Transect, mean))
P_minus = Pmeans[1:13]
P_plus = Pmeans[14:26]
t.test(P_plus, P_minus, paired=TRUE)

    Paired t-test

data:  P_plus and P_minus
t = -0.27532, df = 12, p-value = 0.7878
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.01663756  0.01290460
sample estimates:
mean of the differences 
           -0.001866481 

A linear model is also constructed as both a supplement and a comparison to the paired t-test. We get practically the same p-values, but we also get to see diagnostic plots to confirm that assumptions are met reasonably well.

The linear model is constructed with P (or log(P)) as the dependent variable, transect as the first factor, and microsite as the second factor. We first use the raw P values, not transformed. The residuals look fairly normally distributed, but the first residual plot does not look great.
No evidence for association between P and microsite (p=0.77).

fit <- lm(P ~ Transect + Microsite, data=tissuenutrients)
anova(fit)
Analysis of Variance Table

Response: P
          Df    Sum Sq    Mean Sq F value Pr(>F)
Transect  12 0.0033813 2.8178e-04  1.1784 0.3762
Microsite  1 0.0000212 2.1228e-05  0.0888 0.7698
Residuals 15 0.0035869 2.3912e-04               
drop1(fit, test="F")
Single term deletions

Model:
P ~ Transect + Microsite
          Df Sum of Sq       RSS     AIC F value Pr(>F)
<none>                 0.0035869 -232.94               
Transect  12 0.0034023 0.0069892 -237.59  1.1857 0.3720
Microsite  1 0.0000212 0.0036081 -234.76  0.0888 0.7698
layout(matrix(1:4,1,4, byrow=T))
plot(fit)

We can re-do this analysis, using log(P) values this time. Still no evidence of an association (p=0.71). The residual plots look better (more homogeneous variance), so this analysis is preferred over the analysis of P (untransformed).

fit <- lm(log(P) ~ Transect + Microsite, data=tissuenutrients)
anova(fit)
Analysis of Variance Table

Response: log(P)
          Df  Sum Sq  Mean Sq F value Pr(>F)
Transect  12 0.39191 0.032659  1.1548 0.3901
Microsite  1 0.00406 0.004060  0.1435 0.7101
Residuals 15 0.42421 0.028281               
drop1(fit, test="F")
Single term deletions

Model:
log(P) ~ Transect + Microsite
          Df Sum of Sq     RSS     AIC F value Pr(>F)
<none>                 0.42421 -94.520               
Transect  12   0.39590 0.82011 -99.403  1.1666 0.3831
Microsite  1   0.00406 0.42827 -96.244  0.1435 0.7101
plot(fit)

Tissue Potassium

Tissue Potassium would be expected to increase at clay-rich microsites due to the decreased effect of nutrient leeching presumably resulting in greater K uptake.

data visualization

layout(matrix(1:2,1,2))
hist(tissuenutrients$K, main="", col="tan", breaks=20)
qqnorm(tissuenutrients$K, main = '', ylab="tissue K")
qqline(tissuenutrients$K)

The histogram for tissue potassium concentrations is slightly right-skewed. Q-Q plots were used to test for normality because the Shapiro-Wilk test is not powerful at small sample size (n=31). Sample quantiles do not show noticible deviation from Q-Q plot does not suggest any significant deviations from normality. Four high outliers (large K values) can be seen in the Q-Q plots. they are on all on different transects and represent the two different microsites equally.

subset(tissuenutrients, K>0.7 )
   Microsite Transect   TN        P    K
3          +        2 0.63 0.080000 0.81
6          -        3 0.81 0.097890 0.89
10         -        5 0.55 0.098468 0.90
28         +       13 0.74 0.065631 0.77

Let’s now just visualize the association between potassium and microsite, and also see the pairing by transect (transects 2 and 13 are in red, transects 3 and 5 are in blue).

transect.color = rep("black",13)
transect.color[c(2,13,3, 5)] = c("red", "red", "blue", "blue")
layout(matrix(1:2,1,2))
with(tissuenutrients, interaction.plot(Microsite, Transect, K, legend=FALSE, col=transect.color))
plot(K~Microsite, data=tissuenutrients)

The outliers are much “less” outliers when we consider the log(potassium) values. The Q-Q plot also looks better (though still has a weird high tail):

layout(matrix(1:4,2,2, byrow=T))
with(tissuenutrients, hist(log(K), main="", col="tan", breaks=20))
with(tissuenutrients, qqnorm(log(K), main = '', ylab="tissue log(K)"))
with(tissuenutrients, interaction.plot(Microsite, Transect, log(K), legend=FALSE, col=transect.color))
plot(log(K)~Microsite, data=tissuenutrients)

From this, 2 transects (those in red) have noticibly K higher in the pitcher-present area, and 2 transects (those in blue) have K higher in the pitcher-absent area.

statistical analysis

If we analyze the log values, it seems that (at first sight) we can use a t-test or a 2-way ANOVA-based method. We will double-check the normality of residuals after we get these residuals.

A new vector is created specifically for the paired t-test.

table(tissuenutrients[1:2])
         Transect
Microsite 1 2 3 4 5 6 7 8 9 10 11 12 13
        - 1 1 1 1 1 1 1 1 1  1  4  1  1
        + 1 1 1 1 1 1 1 1 1  1  1  1  1
logKmeans = with(tissuenutrients, tapply(log(K), Microsite:Transect, mean))
K_minus = logKmeans[1:13]
K_plus = logKmeans[14:26]
K_minus  # log values, in fact
       -:1        -:2        -:3        -:4        -:5        -:6 
-0.5276327 -0.4942963 -0.1165338 -0.4620355 -0.1053605 -0.6161861 
       -:7        -:8        -:9       -:10       -:11       -:12 
-0.4462871 -0.4780358 -0.6539265 -0.5447272 -0.6598031 -0.7985077 
      -:13 
-0.7339692 
K_plus
       +:1        +:2        +:3        +:4        +:5        +:6 
-0.4004776 -0.2107210 -0.6539265 -0.5108256 -0.4942963 -0.6348783 
       +:7        +:8        +:9       +:10       +:11       +:12 
-0.4462871 -0.5276327 -0.6348783 -0.6161861 -0.5621189 -0.8675006 
      +:13 
-0.2613648 

A paired t-test is performed on pitcher plant-absent and pitcher plant-present K tissue concentrations.

t.test(K_plus, K_minus, paired = TRUE)

    Paired t-test

data:  K_plus and K_minus
t = -0.19964, df = 12, p-value = 0.8451
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.1684364  0.1401607
sample estimates:
mean of the differences 
            -0.01413787 

K concentrations are not significantly different between pitcher-present and pitcher-absent microsites (p=0.85).
We can do this again with the non-transformed K values, but there is still no evidence of association (p=0.79).

Kmeans = with(tissuenutrients, tapply(K, Microsite:Transect, mean))
K_minus = Kmeans[1:13]
K_plus = Kmeans[14:26]
t.test(K_plus, K_minus, paired=TRUE)

    Paired t-test

data:  K_plus and K_minus
t = -0.26949, df = 12, p-value = 0.7921
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.11705568  0.09128645
sample estimates:
mean of the differences 
            -0.01288462 

A linear model is also constructed as both a supplement and a comparison to the paired t-test. We get practically the same p-values, but we also get to see diagnostic plots to confirm that assumptions are met reasonably well.

The linear model is constructed with K (or log(K)) as the dependent variable, transect as the first factor, and microsite as the second factor. We first use the raw K values, not transformed. The residuals look fairly normally distributed, but the first residual plot does not look great.
No evidence for association between K and microsite (p=0.82).

fit <- lm(K ~ Transect + Microsite, data=tissuenutrients)
anova(fit)
Analysis of Variance Table

Response: K
          Df  Sum Sq  Mean Sq F value Pr(>F)
Transect  12 0.20410 0.017008  1.4033 0.2643
Microsite  1 0.00068 0.000680  0.0561 0.8160
Residuals 15 0.18180 0.012120               
drop1(fit, test="F")
Single term deletions

Model:
K ~ Transect + Microsite
          Df Sum of Sq     RSS     AIC F value Pr(>F)
<none>                 0.18180 -119.09               
Transect  12   0.20460 0.38640 -121.23  1.4068 0.2629
Microsite  1   0.00068 0.18248 -120.98  0.0561 0.8160
layout(matrix(1:4,1,4, byrow=T))
plot(fit)

We can re-do this analysis, using log(K) values this time. Still no evidence of an association (p=0.89). The residual plots look better (more homogeneous variance), so this analysis is preferred over the analysis of K (untransformed).

fit <- lm(log(K) ~ Transect + Microsite, data=tissuenutrients)
anova(fit)
Analysis of Variance Table

Response: log(K)
          Df  Sum Sq  Mean Sq F value Pr(>F)
Transect  12 0.53474 0.044561  1.6578 0.1762
Microsite  1 0.00058 0.000576  0.0214 0.8856
Residuals 15 0.40321 0.026880               
drop1(fit, test="F")
Single term deletions

Model:
log(K) ~ Transect + Microsite
          Df Sum of Sq     RSS     AIC F value Pr(>F)
<none>                 0.40321 -95.992               
Transect  12   0.53394 0.93714 -95.534  1.6553 0.1769
Microsite  1   0.00058 0.40378 -97.951  0.0214 0.8856
plot(fit)

NPK Limitation Diagram

N:P, N:K, and P:K ratios from tissue nutrient measurements are compared to critical ratios put forth by (Olde Venterink et al. 2003) to determine whether plants are P or P + N limited (N:P > 14.5, K:P >3.4), K or K + N limited (N:K > 3.1, K:P < 3.4), or N limited (N:P < 14.5, N:K < 2.1). The relationship between N:P:K ratios and type of nutrient limitation is illustrated in the subsequent triaxial graph.

library(ggtern)
plot <- ggtern(data = tissuenutrients, aes(x = TN, y = P, z = K)) +
  geom_Tisoprop(value=.3226, color='darkgreen') + # T=Top apex: K/(K+N)=.3226 i.e N/K=1/.3226-1=2.1
  geom_Tisoprop(value=.2439, color='lightgreen') + # Top apex: K/(K+N)=.2439 i.e N/K=1/.2439-1=3.1
  geom_Risoprop(value=.93548, color='darkblue') + # Right apex P/(P+N)=0.145: N/P=1/.06452-1=14.5
  geom_Lisoprop(value=.7727273, color='darkred') + # Left apex K/(K+P)=: P/K=1/0.7727273-1=, K/P=1/0.2941=3.4
  geom_point(aes(fill = Microsite),
             size = 2,
             shape = 21,
             color = "black") +
  ggtitle("NPK limitation diagram") +
  labs(fill = "Microsite") +
  theme_rgbw() +
  theme(legend.position = c(0,1),
        legend.justification = c(1, 1))
plot

Critical ratios are represented by vertical lines.

Plants in both pitcher-present and pitcher-absent microsites are most strongly N-limited (?). CA: I would say that they are most strongly P-limited.

About N limitation:

  • All sites have N:P higher than 14.5 (blue isoproportion line). (so none are N limited, correct?).
  • None of the N:K values are higher than 3.1 (light green line), but some are higher than 2.1 (darkgreen line). (these would definitely not be N limited, correct?)

About P limitation: all points are super far from the top apex (100% P), close to the horizontal line across that apex, where there is 0% P. More specifically:

  • all sites have K:P above 3.4 (red line), and
  • all sites have N:P above 4.5 (blue line).

The points are overlapping a lot around the line N:K = 2.1, so it’s hard to see if there are mostly pitcher-present or pitcher-absent sites above/below this line. But let’s tabulate: we get that there are only 4 sites with TN:K above 2.1. Of these 4 sites, 3 are in a pitcher-absent patch (2 along the same transect 11, the third in transect 13), and the 4th is in a pitcher-present area. It goes in the expected direction, but 3 versus 1 is not statistically different from a half/half ratio.

with(tissuenutrients, table(TN/K<2.1, Microsite))
       Microsite
         -  +
  FALSE  3  1
  TRUE  13 12
subset(tissuenutrients, TN/K>2.1)
   Microsite Transect   TN        P    K
1          +        1 1.55 0.140000 0.67
24         -       11 1.12 0.075815 0.50
25         -       11 1.18 0.075889 0.55
29         -       13 1.17 0.098511 0.48
transect.color = rep("black",13)
transect.color[c(2,13,11)] = c("red", "blue","purple")
layout(matrix(1:2,1,2))
with(tissuenutrients, interaction.plot(Microsite, Transect, TN/K, legend=F, col=transect.color))
plot(TN/K ~ Microsite, data=tissuenutrients)

From this we see that in 10 of the 13 transect, TN:K is higher in the pitcher-present area than in the pitcher-absent patch. Does that run counter to expectations? (in the plot above and those below, transect 2 is in red, transect 13 is in blue, and transect 11 is in purple.)

If Phosphorous is limiting, we could analyze the ratios P/xxx. For example:

layout(matrix(1:4,2,2))
with(tissuenutrients, interaction.plot(Microsite, Transect, P/K, legend=F, col=transect.color))
with(tissuenutrients, interaction.plot(Microsite, Transect, P/TN, legend=F, col=transect.color))
with(tissuenutrients, interaction.plot(Microsite, Transect, P/(K+TN), legend=F, col=transect.color))
with(tissuenutrients, interaction.plot(Microsite, Transect, P/(P+K+TN), legend=F, col=transect.color))

We if look at the proportion of P, i.e. P/(P+K+TN), we still don’t find any evidence of association with microsite (pitcher plant present/absent).

fit = lm(P/(P+K+TN) ~ Transect + Microsite, data=tissuenutrients)
drop1(fit, test="F")
Single term deletions

Model:
P/(P + K + TN) ~ Transect + Microsite
          Df  Sum of Sq        RSS     AIC F value  Pr(>F)  
<none>                  0.00095217 -271.40                  
Transect  12 0.00156373 0.00251590 -267.22  2.0528 0.09461 .
Microsite  1 0.00019328 0.00114545 -268.04  3.0448 0.10144  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
layout(matrix(1:4,1,4))
plot(fit)

How about N+P limitation?

layout(matrix(1:2,1,2))
with(tissuenutrients, interaction.plot(Microsite, Transect, (TN+P)/K, legend=T, col=transect.color))
with(tissuenutrients, interaction.plot(Microsite, Transect, (TN+P)/(TN+P+K), legend=F, col=transect.color))

fit = lm((TN+P)/(TN+P+K) ~ Transect + Microsite, data=tissuenutrients)
drop1(fit, test="F")
Single term deletions

Model:
(TN + P)/(TN + P + K) ~ Transect + Microsite
          Df Sum of Sq      RSS     AIC F value Pr(>F)
<none>                 0.086167 -140.74               
Transect  12  0.087773 0.173940 -144.37  1.2733 0.3246
Microsite  1  0.005863 0.092030 -140.84  1.0206 0.3284
layout(matrix(1:4,1,4))
plot(fit)

soilNPKphys = read.csv(file = "soilNPKphys.csv", header = TRUE)
soilNPKphys$Transect = factor(soilNPKphys$Transect) # interpreted as factor, not numerical values
str(soilNPKphys)
'data.frame':   31 obs. of  9 variables:
 $ Microsite : Factor w/ 2 levels "-","+": 2 1 2 1 2 1 2 1 2 1 ...
 $ Transect  : Factor w/ 13 levels "1","2","3","4",..: 1 1 2 2 3 3 4 4 5 5 ...
 $ NO3       : num  0.6 5.2 0.9 0.9 3.8 2.8 3.2 1.5 1 1.7 ...
 $ NH4       : num  2.3 18.9 1.6 2.6 3.1 1.9 7.3 2 2.8 2.3 ...
 $ K         : num  128 321 105 130 147 ...
 $ P         : num  0.9 0.3 0.2 0.3 0.3 0.3 0.4 0.5 0.3 0.9 ...
 $ Clay_Depth: num  17.8 13 12.7 10.2 20 ...
 $ NumTowers : int  11 5 4 8 6 0 5 0 5 0 ...
 $ Height    : num  -0.182 -0.682 -1.183 -1.683 -0.682 ...
tissueNPKphys = read.csv(file = "tissueNPKphys.csv", header = TRUE)
tissueNPKphys$Transect = factor(tissueNPKphys$Transect) # interpreted as factor, not numerical values
str(tissueNPKphys)
'data.frame':   29 obs. of  8 variables:
 $ Microsite : Factor w/ 2 levels "-","+": 2 1 2 1 2 1 2 1 2 1 ...
 $ Transect  : Factor w/ 13 levels "1","2","3","4",..: 1 1 2 2 3 3 4 4 5 5 ...
 $ TN        : num  1.55 1.07 0.63 0.7 1.09 0.81 1.01 0.52 1.08 0.55 ...
 $ K         : num  0.14 0.0869 0.08 0.06 0.0868 ...
 $ P         : num  0.67 0.59 0.81 0.61 0.52 0.89 0.6 0.63 0.61 0.9 ...
 $ Clay_Depth: num  17.8 13 12.7 10.2 20 ...
 $ NumTowers : int  11 5 4 8 6 0 5 0 5 0 ...
 $ Height    : num  -0.182 -0.682 -1.183 -1.683 -0.682 ...

Analysis of Soil Clay Composition

Clay Depth

Clay depth was measured by placing a soil augar with a 23 cm chamber 69 cm into the soil and recording the maximum height of the clay layer. Measurements were recorded at 1) a random location every 15 m along each transect and 2) following the same stratified random schema and coincident with soil and tissue nutrient measurements. Clay depth measurements were not taken concurrent with the PRS soil nutrient samples to prevent the introduction of local drainage. The number of measurements per transect varies with transect length. All distances are relative.

In general, increased clay content at a given microsite is expected to decrease soil and tissue N levels through inhibition of nitrification and reduction in root functionality as a result of local anoxia resulting from poor drainage at these sites. Conversely, P levels are expected to increase with clay content

layout(matrix(1:2,1,2))
hist(soilNPKphys$Clay_Depth, main="", col="tan", breaks=20)
qqnorm(soilNPKphys$Clay_Depth, main = '', ylab="Clay Depth")
qqline(soilNPKphys$Clay_Depth)

The histogram for clay depth concentrations is slightly left-skewed. Q-Q plots were used to test for normality because the Shapiro-Wilk test is not powerful at small sample size (n=31). Sample quantiles do not show noticible deviation from Q-Q plot does not suggest any significant deviations from normality (though it has a wierd high tail). No noticible outliers can be seen in the Q-Q plots.

Let’s now just visualize the association between clay depth and microsite, and also see the pairing by transect.

layout(matrix(1:2,1,2))
with(soilNPKphys, interaction.plot(Microsite, Transect, Clay_Depth, legend=FALSE))
plot(Clay_Depth~Microsite, data=soilNPKphys)

The outliers are much “less” outliers when we consider the log(potassium) values:

layout(matrix(1:4,2,2, byrow=T))
with(soilNPKphys, hist(log(Clay_Depth), main="", col="tan", breaks=20))
with(soilNPKphys, qqnorm(log(Clay_Depth), main = '', ylab="soil log(Clay Depth)"))
with(soilNPKphys, interaction.plot(Microsite, Transect, log(Clay_Depth), legend=FALSE))
plot(log(Clay_Depth)~Microsite, data=soilNPKphys)

There are few large changes in clay depth between microsite but clay depth, in general, decreases on pitcher-present microsites.

statistical analysis

If we analyze the log values, it seems that (at first sight) we can use a t-test or a 2-way ANOVA-based method. We will double-check the normality of residuals after we get these residuals.

At transects 11, 12, and 13, multiple data points (for soil nutrients only) were collected for pitcher plant-absent microsites. It is necessary to average these pitcher plant-absent data points for their respective transects to perform a t-test (which requires a single value per transect per microsite). This averaging is not necessary with a 2-way ANOVA-based approach.

A new vector is created specifically for the paired t-test.

table(soilNPKphys[1:2])
         Transect
Microsite 1 2 3 4 5 6 7 8 9 10 11 12 13
        - 1 1 1 1 1 1 1 1 1  1  4  2  2
        + 1 1 1 1 1 1 1 1 1  1  1  1  1
logClaymeans = with(soilNPKphys, tapply(log(Clay_Depth), Microsite:Transect, mean))
Clay_minus = logClaymeans[1:13]
Clay_plus = logClaymeans[14:26]
Clay_minus  # log values, in fact
     -:1      -:2      -:3      -:4      -:5      -:6      -:7      -:8 
2.566295 2.318458 2.541602 3.170211 2.148559 2.723924 3.209431 3.209431 
     -:9     -:10     -:11     -:12     -:13 
3.209431 2.878074 2.962735 2.958735 2.628579 
Clay_plus
     +:1      +:2      +:3      +:4      +:5      +:6      +:7      +:8 
2.878074 2.541602 2.995857 3.101218 2.613923 1.491780 2.349230 2.913166 
     +:9     +:10     +:11     +:12     +:13 
2.379083 1.625311 3.011606 2.541602 2.541602 

A paired t-test is performed on pitcher plant-absent and pitcher plant-present K soil concentrations.

t.test(Clay_plus, Clay_minus, paired = TRUE)

    Paired t-test

data:  Clay_plus and Clay_minus
t = -1.6223, df = 12, p-value = 0.1307
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.63827946  0.09344678
sample estimates:
mean of the differences 
             -0.2724163 

Clay depths are not significantly different between pitcher-present and pitcher-absent microsites (p=0.13).
We can do this again with the non-transformed clay depth values, but there is still no evidence of association (p=0.10).

Claymeans = with(soilNPKphys, tapply(Clay_Depth, Microsite:Transect, mean))
Clay_minus = Claymeans[1:13]
Clay_plus = Claymeans[14:26]
t.test(Clay_plus, Clay_minus, paired=TRUE)

    Paired t-test

data:  Clay_plus and Clay_minus
t = -1.7628, df = 12, p-value = 0.1034
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -8.3417526  0.8805026
sample estimates:
mean of the differences 
              -3.730625 

A linear model is also constructed as both a supplement and a comparison to the paired t-test. We get practically the same p-values, but we also get to see diagnostic plots to confirm that assumptions are met reasonably well.

The linear model is constructed with clay (or log(clay)) as the dependent variable, transect as the first factor, and microsite as the second factor. We first use the raw clay depth values, not transformed. The residuals look fairly normally distributed, but the first residual plot and Q-Q plot do not look great.
No evidence for association between K and microsite (p=0.09).

fit <- lm(Clay_Depth ~ Transect + Microsite, data=soilNPKphys)
anova(fit)
Analysis of Variance Table

Response: Clay_Depth
          Df Sum Sq Mean Sq F value  Pr(>F)  
Transect  12 460.33  38.361  1.3399 0.28308  
Microsite  1  91.65  91.653  3.2012 0.09141 .
Residuals 17 486.72  28.631                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
drop1(fit, test="F")
Single term deletions

Model:
Clay_Depth ~ Transect + Microsite
          Df Sum of Sq    RSS    AIC F value  Pr(>F)  
<none>                 486.72 113.36                  
Transect  12    430.55 917.27 109.01  1.2532 0.32666  
Microsite  1     91.65 578.37 116.71  3.2012 0.09141 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
layout(matrix(1:4,1,4, byrow=T))
plot(fit)

We can re-do this analysis, using log(clay) values this time. Still no evidence of an association (p=0.11). The residual and Q-Q plots look better (more homogeneous variance), so this analysis is preferred over the analysis of clay (untransformed).

fit <- lm(log(Clay_Depth) ~ Transect + Microsite, data=soilNPKphys)
anova(fit)
Analysis of Variance Table

Response: log(Clay_Depth)
          Df  Sum Sq Mean Sq F value Pr(>F)
Transect  12 2.57569 0.21464  1.3251 0.2901
Microsite  1 0.47465 0.47465  2.9303 0.1051
Residuals 17 2.75368 0.16198               
drop1(fit, test="F")
Single term deletions

Model:
log(Clay_Depth) ~ Transect + Microsite
          Df Sum of Sq    RSS     AIC F value Pr(>F)
<none>                 2.7537 -47.053               
Transect  12   2.38829 5.1420 -51.693  1.2287 0.3400
Microsite  1   0.47465 3.2283 -44.123  2.9303 0.1051
plot(fit)

Next, we evaluate the effect of clay depth on soil nutrients using a linear model.

The linear model is constructed with (or log()) as the dependent variable, transect as the first factor, and clay depth as the second factor. We first use the raw nutrient values, not transformed.

No evidence for association between soil TN and clay depth (p=0.57), soil P and clay depth (p=0.74), or soil K and clay depth (p=0.96).

We can re-do this analysis, using log() values this time. Still no evidence of an association (TN, p=0.85; P, p=0.64; K, p=0.92). The residual and Q-Q plots look better (more homogeneous variance), so this analysis is preferred over the analysis of clay (untransformed).

soilNPKphys$TN <- soilNPKphys$NO3 + soilNPKphys$NH4
fit <- lm(TN ~ Transect + Clay_Depth, data=soilNPKphys)
anova(fit)
Analysis of Variance Table

Response: TN
           Df Sum Sq Mean Sq F value Pr(>F)
Transect   12 244.56  20.380  0.6341 0.7866
Clay_Depth  1  10.68  10.675  0.3321 0.5720
Residuals  17 546.39  32.141               
drop1(fit, test="F")
Single term deletions

Model:
TN ~ Transect + Clay_Depth
           Df Sum of Sq    RSS    AIC F value Pr(>F)
<none>                  546.39 116.95               
Transect   12   254.603 800.99 104.81  0.6601 0.7652
Clay_Depth  1    10.675 557.07 115.55  0.3321 0.5720
plot(fit)

fit <- lm(P ~ Transect + Clay_Depth, data=soilNPKphys)
anova(fit)
Analysis of Variance Table

Response: P
           Df  Sum Sq  Mean Sq F value Pr(>F)
Transect   12 2.17141 0.180951  1.3032 0.3008
Clay_Depth  1 0.01585 0.015852  0.1142 0.7396
Residuals  17 2.36048 0.138852               
drop1(fit, test="F")
Single term deletions

Model:
P ~ Transect + Clay_Depth
           Df Sum of Sq    RSS     AIC F value Pr(>F)
<none>                  2.3605 -51.829               
Transect   12   2.18048 4.5410 -55.546  1.3086 0.2981
Clay_Depth  1   0.01585 2.3763 -53.621  0.1142 0.7396
plot(fit)

fit <- lm(K ~ Transect + Clay_Depth, data=soilNPKphys)
anova(fit)
Analysis of Variance Table

Response: K
           Df Sum Sq Mean Sq F value Pr(>F)
Transect   12  39671  3306.0  0.9998 0.4878
Clay_Depth  1      7     7.2  0.0022 0.9634
Residuals  17  56215  3306.8               
drop1(fit, test="F")
Single term deletions

Model:
K ~ Transect + Clay_Depth
           Df Sum of Sq   RSS    AIC F value Pr(>F)
<none>                  56215 260.59               
Transect   12     39664 95878 253.14  0.9996 0.4880
Clay_Depth  1         7 56222 258.60  0.0022 0.9634
plot(fit)

fit <- lm(log(TN) ~ Transect + Clay_Depth, data=soilNPKphys)
anova(fit)
Analysis of Variance Table

Response: log(TN)
           Df Sum Sq Mean Sq F value Pr(>F)
Transect   12 5.0649 0.42207  0.8801 0.5806
Clay_Depth  1 0.0171 0.01709  0.0356 0.8525
Residuals  17 8.1528 0.47958               
drop1(fit, test="F")
Single term deletions

Model:
log(TN) ~ Transect + Clay_Depth
           Df Sum of Sq     RSS     AIC F value Pr(>F)
<none>                   8.1528 -13.405               
Transect   12    4.7758 12.9286 -23.111  0.8299 0.6220
Clay_Depth  1    0.0171  8.1699 -15.340  0.0356 0.8525
plot(fit)

fit <- lm(log(P) ~ Transect + Clay_Depth, data=soilNPKphys)
anova(fit)
Analysis of Variance Table

Response: log(P)
           Df Sum Sq Mean Sq F value Pr(>F)
Transect   12 5.8384 0.48653  1.6629 0.1643
Clay_Depth  1 0.0656 0.06564  0.2243 0.6418
Residuals  17 4.9740 0.29259               
drop1(fit, test="F")
Single term deletions

Model:
log(P) ~ Transect + Clay_Depth
           Df Sum of Sq     RSS     AIC F value Pr(>F)
<none>                   4.9740 -28.723               
Transect   12    5.8773 10.8513 -28.541  1.6739 0.1612
Clay_Depth  1    0.0656  5.0396 -30.316  0.2243 0.6418
plot(fit)

fit <- lm(log(K) ~ Transect + Clay_Depth, data=soilNPKphys)
anova(fit)
Analysis of Variance Table

Response: log(K)
           Df Sum Sq  Mean Sq F value Pr(>F)
Transect   12 2.1201 0.176677  0.9447 0.5294
Clay_Depth  1 0.0020 0.001986  0.0106 0.9191
Residuals  17 3.1795 0.187027               
drop1(fit, test="F")
Single term deletions

Model:
log(K) ~ Transect + Clay_Depth
           Df Sum of Sq    RSS     AIC F value Pr(>F)
<none>                  3.1795 -42.596               
Transect   12   2.12152 5.3010 -50.749  0.9453 0.5289
Clay_Depth  1   0.00199 3.1814 -44.576  0.0106 0.9191
plot(fit)

We next explore the interaction between tissue nutrient concentrations and clay content using the same linear model design:

The linear model is constructed with clay (or log(clay)) as the dependent variable, transect as the first factor, and microsite as the second factor. We first use the raw clay depth values, not transformed. No evidence for association between tissue TN and clay depth (p=0.57), tissue P and clay depth (p=0.40), or tissue K and clay depth (p=0.36).

We can re-do this analysis, using log() values this time. Still no evidence of an association (TN, p=0.60; P, p=0.44; K, p=0.38). The residual and Q-Q plots look better (more homogeneous variance), so this analysis is preferred over the analysis of clay (untransformed).

fit <- lm(TN ~ Transect + Clay_Depth, data=tissueNPKphys)
anova(fit)
Analysis of Variance Table

Response: TN
           Df  Sum Sq  Mean Sq F value Pr(>F)
Transect   12 0.73255 0.061046  1.4956 0.2282
Clay_Depth  1 0.01365 0.013651  0.3344 0.5716
Residuals  15 0.61227 0.040818               
drop1(fit, test="F")
Single term deletions

Model:
TN ~ Transect + Clay_Depth
           Df Sum of Sq     RSS     AIC F value Pr(>F)
<none>                  0.61227 -83.878               
Transect   12   0.73140 1.34367 -85.085  1.4932 0.2291
Clay_Depth  1   0.01365 0.62592 -85.239  0.3344 0.5716
plot(fit)

fit <- lm(P ~ Transect + Clay_Depth, data=tissueNPKphys)
anova(fit)
Analysis of Variance Table

Response: P
           Df   Sum Sq   Mean Sq F value Pr(>F)
Transect   12 0.204099 0.0170083  1.4681 0.2384
Clay_Depth  1 0.008705 0.0087051  0.7514 0.3997
Residuals  15 0.173775 0.0115850               
drop1(fit, test="F")
Single term deletions

Model:
P ~ Transect + Clay_Depth
           Df Sum of Sq     RSS     AIC F value Pr(>F)
<none>                  0.17377 -120.40               
Transect   12  0.188142 0.36192 -123.12  1.3533 0.2861
Clay_Depth  1  0.008705 0.18248 -120.98  0.7514 0.3997
plot(fit)

fit <- lm(K ~ Transect + Clay_Depth, data=tissueNPKphys)
anova(fit)
Analysis of Variance Table

Response: K
           Df    Sum Sq    Mean Sq F value Pr(>F)
Transect   12 0.0033813 0.00028178  1.2398 0.3421
Clay_Depth  1 0.0001989 0.00019888  0.8750 0.3644
Residuals  15 0.0034092 0.00022728               
drop1(fit, test="F")
Single term deletions

Model:
K ~ Transect + Clay_Depth
           Df Sum of Sq       RSS     AIC F value Pr(>F)
<none>                  0.0034092 -234.41               
Transect   12 0.0034610 0.0068702 -238.09   1.269 0.3268
Clay_Depth  1 0.0001989 0.0036081 -234.76   0.875 0.3644
plot(fit)

fit <- lm(log(TN) ~ Transect + Clay_Depth, data=tissueNPKphys)
anova(fit)
Analysis of Variance Table

Response: log(TN)
           Df  Sum Sq  Mean Sq F value Pr(>F)
Transect   12 0.91612 0.076343  1.4766 0.2352
Clay_Depth  1 0.01463 0.014632  0.2830 0.6025
Residuals  15 0.77555 0.051703               
drop1(fit, test="F")
Single term deletions

Model:
log(TN) ~ Transect + Clay_Depth
           Df Sum of Sq     RSS     AIC F value Pr(>F)
<none>                  0.77555 -77.023               
Transect   12   0.91064 1.68619 -78.500  1.4677 0.2386
Clay_Depth  1   0.01463 0.79018 -78.481  0.2830 0.6025
plot(fit)

fit <- lm(log(P) ~ Transect + Clay_Depth, data=tissueNPKphys)
anova(fit)
Analysis of Variance Table

Response: log(P)
           Df  Sum Sq  Mean Sq F value Pr(>F)
Transect   12 0.53474 0.044561  1.7235 0.1587
Clay_Depth  1 0.01597 0.015966  0.6175 0.4442
Residuals  15 0.38782 0.025855               
drop1(fit, test="F")
Single term deletions

Model:
log(P) ~ Transect + Clay_Depth
           Df Sum of Sq     RSS     AIC F value Pr(>F)
<none>                  0.38782 -97.121               
Transect   12   0.50210 0.88992 -97.034  1.6184 0.1876
Clay_Depth  1   0.01597 0.40378 -97.951  0.6175 0.4442
plot(fit)

fit <- lm(log(K) ~ Transect + Clay_Depth, data=tissueNPKphys)
anova(fit)
Analysis of Variance Table

Response: log(K)
           Df  Sum Sq  Mean Sq F value Pr(>F)
Transect   12 0.39191 0.032659  1.2061 0.3605
Clay_Depth  1 0.02208 0.022078  0.8153 0.3808
Residuals  15 0.40619 0.027079               
drop1(fit, test="F")
Single term deletions

Model:
log(K) ~ Transect + Clay_Depth
           Df Sum of Sq     RSS      AIC F value Pr(>F)
<none>                  0.40619  -95.779               
Transect   12   0.39659 0.80278 -100.022  1.2205 0.3525
Clay_Depth  1   0.02208 0.42827  -96.244  0.8153 0.3808
plot(fit)

Crayfish Burrow Exits as a proxy for clay content

Crayfish, Fallicambarus spp., abundant in wet pine savannas create distinctive soil changes by constructing “chimney” and “mound” structures as antechambers extensive underground burrows. Both structures (refered to collectively as towers), preferentially constructed based on time of year and weather, consist of silt or clay subsoil and therefore serve as a reasonable approximation of clay content (Brewer 1999).

Crayfish tower abundance was recorded coincident with clay depth by counting the number of towers visible within a 1 X 1 m PVC square directly to the right of the transect.

We explore tower abundance as a proxy and a possibly more accurate measure for soil clay content compared to relative clay depth by applying on a Poisson regression using number of towers as the dependent variable and Transect and Microsite as the first and second factors respectively.

library(MASS)
fit <- glm(soilNPKphys$NumTowers ~ soilNPKphys$Transect + soilNPKphys$Microsite, family=poisson)
summary(fit)

Call:
glm(formula = soilNPKphys$NumTowers ~ soilNPKphys$Transect + 
    soilNPKphys$Microsite, family = poisson)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-2.29640  -1.60832  -0.00022   0.86316   1.79748  

Coefficients:
                         Estimate Std. Error z value Pr(>|z|)    
(Intercept)               1.66268    0.27875   5.965 2.45e-09 ***
soilNPKphys$Transect2    -0.28768    0.38188  -0.753 0.451253    
soilNPKphys$Transect3    -0.98083    0.47871  -2.049 0.040473 *  
soilNPKphys$Transect4    -1.16315    0.51235  -2.270 0.023193 *  
soilNPKphys$Transect5    -1.16315    0.51235  -2.270 0.023193 *  
soilNPKphys$Transect6    -0.98083    0.47871  -2.049 0.040473 *  
soilNPKphys$Transect7    -0.69315    0.43301  -1.601 0.109431    
soilNPKphys$Transect8    -0.98083    0.47871  -2.049 0.040473 *  
soilNPKphys$Transect9    -0.37469    0.39167  -0.957 0.338747    
soilNPKphys$Transect10    0.22314    0.33541   0.665 0.505868    
soilNPKphys$Transect11    0.09524    0.30794   0.309 0.757105    
soilNPKphys$Transect12  -19.23701 1977.86986  -0.010 0.992240    
soilNPKphys$Transect13   -1.95885    0.62989  -3.110 0.001872 ** 
soilNPKphys$Microsite+    0.71004    0.18392   3.861 0.000113 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 132.713  on 30  degrees of freedom
Residual deviance:  49.457  on 17  degrees of freedom
AIC: 152.51

Number of Fisher Scoring iterations: 15

Number of towers shows significant interaction with microsite.

Height as a possible controller of nutrient limitation

It is possible that since elevated microsites would be subject to greater leeching and would therefore be more nutrient poor compared to microsites lower in the landscape irrespective of soil composition. Height relative to the beginning of the transect was approximated with a range finder according to standard practices. Measurements were taken following the same stratified random schema and coincident as soil and tissue nutrient measurements.

layout(matrix(1:2,1,2))
hist(soilNPKphys$Height, main="", col="tan", breaks=20)
qqnorm(soilNPKphys$Height, main = '', ylab="Height")
qqline(soilNPKphys$Height)

The histogram for height concentrations is very slightly right-skewed. Q-Q plots were used to test for normality because the Shapiro-Wilk test is not powerful at small sample size (n=31). Sample quantiles do not show noticible deviation from Q-Q plot does not suggest any significant deviations from normality. No noticible outliers can be seen in the Q-Q plots.

Let’s now just visualize the association between height and microsite, and also see the pairing by transect.

layout(matrix(1:2,1,2))
with(soilNPKphys, interaction.plot(Microsite, Transect, Height, legend=FALSE))
plot(Height~Microsite, data=soilNPKphys)

statistical analysis

If we analyze the log values, it seems that (at first sight) we can use a t-test or a 2-way ANOVA-based method. We will double-check the normality of residuals after we get these residuals.

At transects 11, 12, and 13, multiple data points (for soil nutrients only) were collected for pitcher plant-absent microsites. It is necessary to average these pitcher plant-absent data points for their respective transects to perform a t-test (which requires a single value per transect per microsite). This averaging is not necessary with a 2-way ANOVA-based approach.

A new vector is created specifically for the paired t-test.

Height is not significantly different between pitcher-present and pitcher-absent microsites (p=0.11).

Heightmeans = with(soilNPKphys, tapply(Height, Microsite:Transect, mean))
Height_minus = Heightmeans[1:13]
Height_plus = Heightmeans[14:26]
t.test(Height_plus, Height_minus, paired=TRUE)

    Paired t-test

data:  Height_plus and Height_minus
t = 1.7339, df = 12, p-value = 0.1085
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.07649503  0.67264888
sample estimates:
mean of the differences 
              0.2980769 

A linear model is also constructed as both a supplement and a comparison to the paired t-test. We get practically the same p-values, but we also get to see diagnostic plots to confirm that assumptions are met reasonably well.

The linear model is constructed with height as the dependent variable, transect as the first factor, and microsite as the second factor. We first use the raw height values, not transformed.

No evidence for association between K and microsite (p=0.24).

fit <- lm(Height ~ Transect + Microsite, data=soilNPKphys)
anova(fit)
Analysis of Variance Table

Response: Height
          Df Sum Sq Mean Sq F value    Pr(>F)    
Transect  12 36.929 3.07740  7.6375 0.0001082 ***
Microsite  1  0.608 0.60845  1.5100 0.2358714    
Residuals 17  6.850 0.40293                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
drop1(fit, test="F")
Single term deletions

Model:
Height ~ Transect + Microsite
          Df Sum of Sq    RSS     AIC F value   Pr(>F)    
<none>                  6.850 -18.802                     
Transect  12    37.322 44.172  14.977  7.7188 0.000101 ***
Microsite  1     0.608  7.458 -18.164  1.5100 0.235871    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
layout(matrix(1:4,1,4, byrow=T))
plot(fit)

Next, we evaluate the effect of height on soil nutrients using a linear model.

The linear model is constructed with (or log()) as the dependent variable, transect as the first factor, and height as the second factor.

We first use the raw nutrient values, not transformed.

No evidence for association between soil TN and height (p=0.57), soil P and height (p=0.74), or soil K and height (p=0.96).

We can re-do this analysis, using log() values this time. Still no evidence of an association (TN, p=0.85; P, p=0.64; K, p=0.92). The residual and Q-Q plots look better (more homogeneous variance), so this analysis is preferred over the analysis of height (untransformed).

fit <- lm(TN ~ Transect + Height, data=soilNPKphys)
anova(fit)
Analysis of Variance Table

Response: TN
          Df Sum Sq Mean Sq F value Pr(>F)
Transect  12 244.56  20.380  0.6489 0.7744
Height     1  23.19  23.185  0.7383 0.4022
Residuals 17 533.88  31.405               
drop1(fit, test="F")
Single term deletions

Model:
TN ~ Transect + Height
         Df Sum of Sq    RSS    AIC F value Pr(>F)
<none>                533.88 116.23               
Transect 12   264.680 798.56 104.71  0.7023 0.7298
Height    1    23.185 557.07 115.55  0.7383 0.4022
plot(fit)

fit <- lm(P ~ Transect + Height, data=soilNPKphys)
anova(fit)
Analysis of Variance Table

Response: P
          Df  Sum Sq  Mean Sq F value Pr(>F)
Transect  12 2.17141 0.180951  1.3068 0.2990
Height     1 0.02236 0.022356  0.1614 0.6928
Residuals 17 2.35398 0.138469               
drop1(fit, test="F")
Single term deletions

Model:
P ~ Transect + Height
         Df Sum of Sq    RSS     AIC F value Pr(>F)
<none>                2.3540 -51.914               
Transect 12   1.72934 4.0833 -58.839  1.0407 0.4583
Height    1   0.02236 2.3763 -53.621  0.1614 0.6928
plot(fit)

fit <- lm(K ~ Transect + Height, data=soilNPKphys)
anova(fit)
Analysis of Variance Table

Response: K
          Df Sum Sq Mean Sq F value  Pr(>F)  
Transect  12  39671  3306.0  1.2092 0.35100  
Height     1   9742  9742.2  3.5632 0.07627 .
Residuals 17  46480  2734.1                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
drop1(fit, test="F")
Single term deletions

Model:
K ~ Transect + Height
         Df Sum of Sq   RSS    AIC F value  Pr(>F)  
<none>                46480 254.70                  
Transect 12     43201 89681 251.07  1.3167 0.29415  
Height    1      9742 56222 258.60  3.5632 0.07627 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(fit)

fit <- lm(log(TN) ~ Transect + Height, data=soilNPKphys)
anova(fit)
Analysis of Variance Table

Response: log(TN)
          Df Sum Sq Mean Sq F value Pr(>F)
Transect  12 5.0649 0.42207  0.9013 0.5635
Height     1 0.2093 0.20930  0.4470 0.5128
Residuals 17 7.9606 0.46827               
drop1(fit, test="F")
Single term deletions

Model:
log(TN) ~ Transect + Height
         Df Sum of Sq     RSS     AIC F value Pr(>F)
<none>                 7.9606 -14.144               
Transect 12    5.2281 13.1886 -22.494  0.9304 0.5405
Height    1    0.2093  8.1699 -15.340  0.4470 0.5128
plot(fit)

fit <- lm(log(P) ~ Transect + Height, data=soilNPKphys)
anova(fit)
Analysis of Variance Table

Response: log(P)
          Df Sum Sq Mean Sq F value Pr(>F)
Transect  12 5.8384 0.48653  1.6458 0.1691
Height     1 0.0141 0.01410  0.0477 0.8297
Residuals 17 5.0256 0.29562               
drop1(fit, test="F")
Single term deletions

Model:
log(P) ~ Transect + Height
         Df Sum of Sq    RSS     AIC F value Pr(>F)
<none>                5.0256 -28.403               
Transect 12    4.8323 9.8578 -31.517  1.3622 0.2727
Height    1    0.0141 5.0396 -30.316  0.0477 0.8297
plot(fit)

fit <- lm(log(K) ~ Transect + Height, data=soilNPKphys)
anova(fit)
Analysis of Variance Table

Response: log(K)
          Df  Sum Sq Mean Sq F value Pr(>F)
Transect  12 2.12012 0.17668  1.0747 0.4349
Height     1 0.38671 0.38671  2.3523 0.1435
Residuals 17 2.79473 0.16440               
drop1(fit, test="F")
Single term deletions

Model:
log(K) ~ Transect + Height
         Df Sum of Sq    RSS     AIC F value Pr(>F)
<none>                2.7947 -46.594               
Transect 12   2.12893 4.9237 -53.038  1.0792 0.4319
Height    1   0.38671 3.1814 -44.576  2.3523 0.1435
plot(fit)

We next explore the interaction between tissue nutrient concentrations and height using the same linear model design:

The linear model is constructed with ) values this time. Still no evidence of an association (TN, p=0.87; P, p=0.29; K, p=0.53). The residual and Q-Q plots look better (more homogeneous variance), so this analysis is preferred over the analysis of height (untransformed).

fit <- lm(TN ~ Transect + Height, data=tissueNPKphys)
anova(fit)
Analysis of Variance Table

Response: TN
          Df  Sum Sq  Mean Sq F value Pr(>F)
Transect  12 0.73255 0.061046  1.4680 0.2385
Height     1 0.00216 0.002155  0.0518 0.8230
Residuals 15 0.62376 0.041584               
drop1(fit, test="F")
Single term deletions

Model:
TN ~ Transect + Height
         Df Sum of Sq     RSS     AIC F value Pr(>F)
<none>                0.62376 -83.339               
Transect 12   0.73455 1.35831 -84.771  1.4720 0.2369
Height    1   0.00216 0.62592 -85.239  0.0518 0.8230
plot(fit)

fit <- lm(P ~ Transect + Height, data=tissueNPKphys)
anova(fit)
Analysis of Variance Table

Response: P
          Df  Sum Sq   Mean Sq F value Pr(>F)
Transect  12 0.20410 0.0170083  1.4780 0.2347
Height     1 0.00987 0.0098698  0.8577 0.3690
Residuals 15 0.17261 0.0115073               
drop1(fit, test="F")
Single term deletions

Model:
P ~ Transect + Height
         Df Sum of Sq     RSS     AIC F value Pr(>F)
<none>                0.17261 -120.60               
Transect 12   0.20314 0.37575 -122.04  1.4711 0.2373
Height    1   0.00987 0.18248 -120.98  0.8577 0.3690
plot(fit)

fit <- lm(K ~ Transect + Height, data=tissueNPKphys)
anova(fit)
Analysis of Variance Table

Response: K
          Df    Sum Sq    Mean Sq F value Pr(>F)
Transect  12 0.0033813 2.8178e-04  1.1902 0.3694
Height     1 0.0000570 5.6983e-05  0.2407 0.6308
Residuals 15 0.0035511 2.3674e-04               
drop1(fit, test="F")
Single term deletions

Model:
K ~ Transect + Height
         Df Sum of Sq       RSS     AIC F value Pr(>F)
<none>                0.0035511 -233.23               
Transect 12 0.0034226 0.0069737 -237.65  1.2048 0.3612
Height    1 0.0000570 0.0036081 -234.76  0.2407 0.6308
plot(fit)

fit <- lm(log(TN) ~ Transect + Height, data=tissueNPKphys)
anova(fit)
Analysis of Variance Table

Response: log(TN)
          Df  Sum Sq  Mean Sq F value Pr(>F)
Transect  12 0.91612 0.076343  1.4520 0.2446
Height     1 0.00149 0.001488  0.0283 0.8686
Residuals 15 0.78869 0.052580               
drop1(fit, test="F")
Single term deletions

Model:
log(TN) ~ Transect + Height
         Df Sum of Sq     RSS     AIC F value Pr(>F)
<none>                0.78869 -76.536               
Transect 12   0.91695 1.70564 -78.167  1.4533 0.2441
Height    1   0.00149 0.79018 -78.481  0.0283 0.8686
plot(fit)

fit <- lm(log(P) ~ Transect + Height, data=tissueNPKphys)
anova(fit)
Analysis of Variance Table

Response: log(P)
          Df  Sum Sq  Mean Sq F value Pr(>F)
Transect  12 0.53474 0.044561  1.7900 0.1428
Height     1 0.03036 0.030360  1.2195 0.2869
Residuals 15 0.37342 0.024895               
drop1(fit, test="F")
Single term deletions

Model:
log(P) ~ Transect + Height
         Df Sum of Sq     RSS     AIC F value Pr(>F)
<none>                0.37342 -98.218               
Transect 12   0.52666 0.90008 -96.704  1.7630 0.1490
Height    1   0.03036 0.40378 -97.951  1.2195 0.2869
plot(fit)

fit <- lm(log(K) ~ Transect + Height, data=tissueNPKphys)
anova(fit)
Analysis of Variance Table

Response: log(K)
          Df  Sum Sq  Mean Sq F value Pr(>F)
Transect  12 0.39191 0.032659  1.1754 0.3780
Height     1 0.01149 0.011494  0.4137 0.5298
Residuals 15 0.41677 0.027785               
drop1(fit, test="F")
Single term deletions

Model:
log(K) ~ Transect + Height
         Df Sum of Sq     RSS     AIC F value Pr(>F)
<none>                0.41677 -95.033               
Transect 12   0.40141 0.81818 -99.471  1.2039 0.3617
Height    1   0.01149 0.42827 -96.244  0.4137 0.5298
plot(fit)

Species Distribution

Species richness and abundance were measured by estimating percent ground cover of each species along each transect 1) at a random location every 15 m along each transect coincident with clay depth measurements and 2) following the same stratified random schema and coincident soil and tissue nutrient measurements.

Species distribution and productivity is influenced by nutrient limitation, with some plants having special adaptations (e.g. symbiotic N2 fixation, root secretion of phosphatases, etc.) that favor nutrient-poor microsites (Olde Venterink et al. 2003).

Clay content, as a proxy for local anoxia, may also influence species distribution according to species flood tolerance.

We use Nonmetric Multidimensional scaling as the ordination method.

First, we look at the species distributions on sites coincident with clay depth data:

speciesClay = read.csv(file = "speciesAbundanceClay.csv", header = TRUE)
library(vegan)
Loading required package: permute
Loading required package: lattice
This is vegan 2.4-3
ord <- metaMDS(speciesClay[9:35], distance = "jaccard",k=2,trymax=1000,autotransform=TRUE,expand=FALSE, plot=FALSE)
Run 0 stress 0.1720377 
Run 1 stress 0.1923467 
Run 2 stress 0.1730232 
Run 3 stress 0.172013 
... New best solution
... Procrustes: rmse 0.01030351  max resid 0.07533652 
Run 4 stress 0.1988807 
Run 5 stress 0.1717517 
... New best solution
... Procrustes: rmse 0.02555623  max resid 0.1629871 
Run 6 stress 0.1730221 
Run 7 stress 0.1921792 
Run 8 stress 0.1996325 
Run 9 stress 0.1720869 
... Procrustes: rmse 0.02679115  max resid 0.1664139 
Run 10 stress 0.1730208 
Run 11 stress 0.1726898 
Run 12 stress 0.1730249 
Run 13 stress 0.1851641 
Run 14 stress 0.1727169 
Run 15 stress 0.1730269 
Run 16 stress 0.1720881 
... Procrustes: rmse 0.02686078  max resid 0.166446 
Run 17 stress 0.1721017 
... Procrustes: rmse 0.02744538  max resid 0.1667545 
Run 18 stress 0.1729556 
Run 19 stress 0.1724165 
Run 20 stress 0.1908043 
Run 21 stress 0.196815 
Run 22 stress 0.1717445 
... New best solution
... Procrustes: rmse 0.001226186  max resid 0.006088162 
... Similar to previous best
*** Solution reached
scl <- 3 ## scaling == 3
colvec <- c("red2", "mediumblue")
plot(ord, type = "n")
with(speciesClay, points(ord, display = "sites", col = colvec,
                      pch = 21, bg = colvec))
text(ord, display = "sites", cex = 0.8, col = "darkcyan")
with(speciesClay, legend("topleft", levels(Microsite), bty = "n",
                      col = colvec, pch = 21, pt.bg = colvec))

adonis(speciesClay[9:35]~speciesClay$Microsite,method="jaccard")

Call:
adonis(formula = speciesClay[9:35] ~ speciesClay$Microsite, method = "jaccard") 

Permutation: free
Number of permutations: 999

Terms added sequentially (first to last)

                      Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)    
speciesClay$Microsite  1     2.225 2.22499  7.3231 0.11042  0.001 ***
Residuals             59    17.926 0.30383         0.88958           
Total                 60    20.151                 1.00000           
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
adonis(speciesClay[11:31]~speciesClay$Clay_cm,method="jaccard")

Call:
adonis(formula = speciesClay[11:31] ~ speciesClay$Clay_cm, method = "jaccard") 

Permutation: free
Number of permutations: 999

Terms added sequentially (first to last)

                    Df SumsOfSqs MeanSqs F.Model     R2 Pr(>F)
speciesClay$Clay_cm  1    0.2966 0.29656 0.89232 0.0149  0.544
Residuals           59   19.6085 0.33235         0.9851       
Total               60   19.9051                 1.0000       

Species show significant grouping by microsite (p=0.001) but not by clay depth (p=0.793)

Next, we look at the species distributions on sites coincident with clay depth data:

speciesNPK = read.csv(file = "speciesAbundanceNPK.csv", header = TRUE)
library(vegan)
ord <- metaMDS(speciesNPK[11:31], distance = "jaccard",k=2,trymax=1000,autotransform=TRUE,expand=FALSE, plot=FALSE)
Run 0 stress 0.1658997 
Run 1 stress 0.1742643 
Run 2 stress 0.1615547 
... New best solution
... Procrustes: rmse 0.0727493  max resid 0.2859159 
Run 3 stress 0.180824 
Run 4 stress 0.1895995 
Run 5 stress 0.1606823 
... New best solution
... Procrustes: rmse 0.04366153  max resid 0.1825723 
Run 6 stress 0.173678 
Run 7 stress 0.1673051 
Run 8 stress 0.1791267 
Run 9 stress 0.1791192 
Run 10 stress 0.1592366 
... New best solution
... Procrustes: rmse 0.05993994  max resid 0.2076739 
Run 11 stress 0.1592368 
... Procrustes: rmse 0.0002579027  max resid 0.0007659765 
... Similar to previous best
Run 12 stress 0.1697659 
Run 13 stress 0.1654879 
Run 14 stress 0.1649254 
Run 15 stress 0.1615097 
Run 16 stress 0.1637516 
Run 17 stress 0.1673051 
Run 18 stress 0.1607993 
Run 19 stress 0.1685889 
Run 20 stress 0.1688124 
*** Solution reached
scl <- 3 ## scaling == 3
colvec <- c("red2", "mediumblue")
plot(ord, type = "n")
with(speciesNPK, points(ord, display = "sites", col = colvec,
                      pch = 21, bg = colvec))
text(ord, display = "sites", cex = 0.8, col = "darkcyan")
with(speciesNPK, legend("topleft", levels(Microsite), bty = "n",
                      col = colvec, pch = 21, pt.bg = colvec))

adonis(speciesNPK[11:31]~speciesNPK$Microsite,method="jaccard")

Call:
adonis(formula = speciesNPK[11:31] ~ speciesNPK$Microsite, method = "jaccard") 

Permutation: free
Number of permutations: 999

Terms added sequentially (first to last)

                     Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)    
speciesNPK$Microsite  1    0.8521 0.85212  3.9075 0.11874  0.001 ***
Residuals            29    6.3241 0.21807         0.88126           
Total                30    7.1762                 1.00000           
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
adonis(speciesNPK[11:31]~speciesNPK$NH4,method="jaccard")

Call:
adonis(formula = speciesNPK[11:31] ~ speciesNPK$NH4, method = "jaccard") 

Permutation: free
Number of permutations: 999

Terms added sequentially (first to last)

               Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)  
speciesNPK$NH4  1    0.3941 0.39415  1.6854 0.05492  0.068 .
Residuals      29    6.7821 0.23387         0.94508         
Total          30    7.1762                 1.00000         
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
adonis(speciesNPK[11:31]~speciesNPK$Clay_Depth,method="jaccard")

Call:
adonis(formula = speciesNPK[11:31] ~ speciesNPK$Clay_Depth, method = "jaccard") 

Permutation: free
Number of permutations: 999

Terms added sequentially (first to last)

                      Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)
speciesNPK$Clay_Depth  1    0.1646 0.16463 0.68089 0.02294  0.802
Residuals             29    7.0116 0.24178         0.97706       
Total                 30    7.1762                 1.00000       

Species do not show significant grouping by clay depth (p=.786) or ammonium (p=0.071) but do show significant grouping by microsite (p=0.001)

R version 3.3.2 (2016-10-31)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 16.04.1 LTS

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] vegan_2.4-3     lattice_0.20-34 permute_0.9-4   MASS_7.3-47    
[5] ggtern_2.2.0    ggplot2_2.2.1  

loaded via a namespace (and not attached):
 [1] latex2exp_0.4.0     Rcpp_0.12.10        DEoptimR_1.0-8     
 [4] plyr_1.8.4          compositions_1.40-1 tools_3.3.2        
 [7] boot_1.3-18         digest_0.6.12       nlme_3.1-129       
[10] evaluate_0.10       tibble_1.3.3        gtable_0.2.0       
[13] mgcv_1.8-16         rlang_0.1.1         Matrix_1.2-7.1     
[16] yaml_2.1.14         parallel_3.3.2      proto_1.0.0        
[19] gridExtra_2.2.1     cluster_2.0.5       stringr_1.2.0      
[22] knitr_1.15.1        rprojroot_1.2       grid_3.3.2         
[25] robustbase_0.92-7   bayesm_3.0-2        rmarkdown_1.4      
[28] tensorA_0.36        magrittr_1.5        backports_1.0.5    
[31] scales_0.4.1        htmltools_0.3.5     colorspace_1.3-2   
[34] labeling_0.3        stringi_1.1.3       energy_1.7-0       
[37] lazyeval_0.2.0      munsell_0.4.3      

References

Brewer, J Stephen. 1999. “Effects of Fire, Competition and Soil Disturbances on Regeneration of a Carnivorous Plant (Drosera Capillaris).” The American Midland Naturalist 141 (1). BioOne: 28–42.

Craft, Christopher B, and Connie Chiang. 2002. “Forms and Amounts of Soil Nitrogen and Phosphorus Across a Longleaf Pine–depressional Wetland Landscape.” Soil Science Society of America Journal 66 (5). Soil Science Society: 1713–21.

Inglett, PW, KR Reddy, and R Corstanje. 2005. “Anaerobic Soils.” Encyclopedia of Soils in the Environment 1. Elsevier FL.

Olde Venterink, H, MJ Wassen, AWM Verkroost, and PC De Ruiter. 2003. “Species Richness–productivity Patterns Differ Between N-, P-, and K-Limited Wetlands.” Ecology 84 (8). Wiley Online Library: 2191–9.